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ABSTRACT 


Remote- sensing  analysis  is  conducted  for  the  Naval  Postgraduate  School  campus, 
containing  buildings,  impervious  surfaces  (asphalt  and  concrete),  natural  ground,  and 
vegetation.  Data  is  from  the  Optech  Titan,  providing  three-wavelength  laser  data  (532, 
1064,  and  1550  nm)  at  10-15  points/m".  Analysis  techniques  for  laser-scanner  (LiDAR) 
data  traditionally  use  only  x,  y,  z  coordinate  information.  The  traditional  approach  is  used 
to  initialize  the  classification  process  into  broad-spatial  classes  (unclassified,  ground, 
vegetation,  and  buildings). 

Spectral  analysis  contributes  a  unique  approach  to  the  classification  process. 
Tools  and  techniques  designed  for  multispectral  imagery  are  adapted  to  the  LiDAR 
analysis  herein.  ENVI’s  N-Dimensional  Visualizer  is  employed  to  develop  training  sets 
for  supervised  classification  techniques,  primarily  Maximum  Likelihood.  Unsupervised 
classification  for  the  combined  spatial/spectral  data  is  accomplished  using  a  K- means 
classifier  for  comparison. 

The  campus  is  classified  into  10  and  16  classes,  compared  to  the  four  from 
traditional  methods.  Addition  of  the  spectral  component  improves  the  discrimination 
among  impervious  surfaces,  other  ground  elements,  and  building  materials.  Maximum 
Likelihood  demonstrates  75%  overall  classification  accuracy,  with  grass  (99.9%),  turf 
(95%),  asphalt  shingles  (94%),  light-building  concrete  (89%),  sand  (88%),  shrubs  (85%), 
asphalt  (84%),  trees  (80%),  and  clay-tile  shingles  (77%).  Post-process  filtering  by 
“number  of  returns”  increases  overall  accuracy  to  82%. 
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I.  INTRODUCTION 

A.  PURPOSE  OF  RESEARCH 

Remote  sensing  from  the  air  and  space  contributes  invaluable  tactical-, 
operational-,  and  strategic-level  intelligence  for  military  commanders  and  national 
leaders.  This  is  especially  true  for  areas  that  are  becoming  increasingly  difficult  for 
United  States’  land-  and  sea-based  platforms  to  access.  Remote- sensing  technologies 
must  continue  to  evolve  in  order  to  maintain  dominance  over  our  adversaries’  expanding 
methods  of  anti-access  area  denial  (A2AD)  and  deception  tactics. 

This  research  is  devoted  to  the  expanding  use  of  the  laser  for  remote- sensing 
purposes  by  analyzing  a  new  and  different  approach  for  terrain  classification  using  data 
from  an  aerial  laser  scanner.  Modem-day  light  detection  and  ranging  (LiDAR)  systems 
offer  more  than  traditional  laser  detection  and  ranging.  Aerial  LiDAR  platforms  are  now 
capable  of  high  resolution,  three-dimensional  (3D),  spatial  and  spectral  imaging  from 
ever-increasing  altitudes.  In  addition  to  active  nighttime  imaging,  LiDAR  gives  us 
exceptional  tree  canopy  and  vegetation  penetration  at  substantially  better  resolutions  than 
synthetic  aperture  radar  (SAR).  This  helps  negate  adversarial  attempts  to  disguise 
buildings  and  vehicles  under  dense  foliage,  and  it  provides  a  means  to  detect  hidden 
roads  and  trails,  even  in  jungles.  LiDAR  is  the  only  remote- sensing  technology  that  can 
regularly  detect  powerlines  and  cables,  a  necessity  for  safe  helicopter  and  unmanned 
aerial  vehicle  (UAV)  operations  in  hostile  and  unfamiliar  terrains.  It  may  even  be  able  to 
detect  the  thin  trip  wires  used  for  some  improvised  explosive  devices  (IED). 

Currently,  the  majority  of  terrain  classification  within  LiDAR  point  clouds  is 
accomplished  using  spatial  attributes,  such  as  point  position  (the  physical  x,  y,  z 
coordinates),  height  above  the  calculated  ground  level,  and  the  apparent  surface 
roughness.  Intensity  is  rarely  included,  and  typically  it  is  only  available  for  a  single 
wavelength.  This  results  in  a  limited  number  of  generalized  classification  categories, 
approximately  two  to  four.  Additional  research  needs  to  analyze  the  conjoined  spectral 
intensities  of  a  multi-wavelength  LiDAR  system  to  determine  if  the  additional  spectral 
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information  collected  from  this  sensor  lends  itself  to  producing  a  greater  variety  of 
classification  results  over  purely  spatial  and  single-intensity  methods. 

B.  OBJECTIVE 

The  primary  objective  of  this  thesis  is  to  determine  if  a  multi-wavelength 
(spectral)  LiDAR  system,  such  as  the  Optech  Titan,  can  generate  more  accurate  and 
diverse  results  in  the  process  of  terrain,  building,  and  vegetation  classification  within  the 
LiDAR  point  cloud.  The  goal  is  to  expand  upon  the  output  of  a  popular  spatial  classifier 
(LASclassify)  that  only  considers  the  geometric  attributes  of  points,  resulting  in  four 
generalized  classes:  unclassified,  ground,  buildings,  and  vegetation.  Herein,  we  will 
analyze  three- wavelength  intensity  values  for  a  5%  random  subset  of  the  point  cloud  to 
distinguish  subclasses  that  correspond  to  specific  manmade  materials  and  terrain  types, 
such  as  asphalt,  concrete,  clay/brick,  turf  (fake  grass),  dirt,  sand,  and  lawn  grass.  Then, 
we  compare  the  effectiveness  of  the  Maximum  Likelihood  (ML),  Spectral  Angle  Mapper 
(SAM),  and  K-means  (Bayesian)  classifiers  in  classifying  the  remaining  95%  of  the  data. 
To  meet  this  objective  in  a  well-understood  environment,  we  select  Optech  Titan  data 
collected  in  three  wavelengths  simultaneously  over  the  campus  of  the  Naval  Postgraduate 
School  in  Monterey,  California — a  semi-urban  (urban  forest)  environment. 

The  secondary  objective  of  this  thesis  is  to  provide  a  detailed  historical  and 
contemporary  account  of  LiDAR  systems  in  space  and  a  thorough  review  of  the 
accompanying  scientific  literature.  One  will  see  that  remote- sensing  laser  systems  have 
been  an  integral  part  of  space  operations  since  the  mid-1960s.  This  secondary  objective 
culminates  in  the  Appendix  where  three  NASA-released  datasets  are  adapted  to  produce 
intriguing  LiDAR  point  clouds  of  the  Moon  and  Mars.  Also,  the  work  detailed  in  the 
Appendix  predates  the  primary  thesis  effort  and  represents  important  practice  with 
LiDAR  data  in  American  Standard  Code  for  Information  Interchange  (ASCII)  format 
using  two  fundamental  software  programs  for  point  cloud  modeling. 
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II.  HISTORICAL  AND  LITERARY  REVIEW  OF  SPACEBORNE 

LiDAR  SYSTEMS 


A.  BACKGROUND 

The  National  Aeronautics  and  Space  Administration  (NASA)  maintains  a  robust 
resume  of  space-based  LiDAR  that  spans  over  five  decades.  Appropriately,  Section  B 
reviews  the  prominent  LiDAR-wielding  space  platforms  designed,  launched,  and 
operated  by  NASA;  Section  B  also  includes  systems  designed,  launched,  and  operated  by 
other  agencies  that  partnered  with  NASA  during  their  endeavors.  Thirteen  systems  are 
presented  in  a  roughly  chronological  order  from  simple  satellites  bearing  only  laser 
reflectors  to  complex  explorer  probes  still  operational  at  the  time  of  writing. 

Section  C  reviews  two  space-based  LiDAR  systems  designed,  launched,  and 
operated  by  nations  other  than  the  United  States  (U.S.):  the  Russian  Space  Program  and 
the  Japan  Aerospace  Exploration  Agency  (JAXA).  In  Section  D,  this  review  concludes 
with  a  look  at  four  new  technologies  and  their  implications  on  the  future  of  LiDAR  in 
space.  Overall,  the  author  located  22  LiDAR- satellite  platforms,  including  the  14  U.S., 
one  French/Russian,  and  the  one  Japanese  sensors  discussed  herein.  Although  not 
covered  further,  additional  noteworthy  space  LiDAR  systems  are  listed,  with  red 
indicating  an  early  mission  failure: 

U.S.  Mars  Polar  Lander  (Atmospheric  LiDAR  for  ice  and  aerosols) 

Mars  Observer  Laser  Altimeter  (Topographic  LiDAR) 

Russia  BALKAN- 1  (Atmospheric  LiDAR  on  the  MIR  space  station) 

BALKAN-2  (Atmospheric  LiDAR  on  the  ALMAZ- IB  satellite) 

China  Chang ’e  Lunar  Orbiter  (Laser  Altimeter) 

Japan  Selenological  and  Engineering  Explorer/Lunar  Orbiter  “SELENE 

aka  Kaguya”  (Laser  Altimeter) 

Table  1  offers  a  comprehensive  list  of  the  upcoming  laser  systems  and  some 
important  facts  for  a  comparison  on  the  evolution  of  space-based  LiDAR.  The  sources  of 
platform  and  orbital  characteristics  are  cited  in  the  narrative  section  for  each  listed 
LiDAR  system.  The  past,  present,  and  future  trends  are  the  author’s  interpretations  based 
on  the  information  in  Table  1. 
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Table  1.  Summary  of  Space-based  LiDAR  Systems  in  Chapter  II. 


Past  (Systems  no  longer  in  operal 

tion) 

Name 

Purpose 

Operating  X 

Beams 

Type 

Orbit 

Beacon 
Explorer-  B/C 

Ground  laser  tracking 
on  satellite  target 

694  nm  (red) 

1 

Ruby 

Elliptical,  low- 
Earth,  79.7° 

Apollo  15  -  17 

Lunar  laser  altimeter 

Unsure  - 
maybe  694  nm 

1 

Q- switch 
Ruby 

Circular  lunar 

ALISSA 

Cloud-profile  LiDAR 
on  Mir  station 

532  nm 
(green) 

1(4 

combined) 

4  Nd:YAG 

Circular,  low- 
Earth,  52° 

Clementine 

Lunar  laser 
mapping — low  quality 

532,  1064  nm 
(IR) 

1 — Either 
separately 

Cr:Nd:YAG 

Circular  lunar  & 
Earth  transfer 

Shuttle  LITE 

Cloud  profile, 
atmospheric  particles 

355  (UV),  532, 
1064  nm 

1—3  Us 
combined 

2  Nd:YAG 

Shuttle  STS -64, 
near-Earth 

Shuttle  Laser 
Altimetry 

Earth  laser 
topography 

1064  nm 

1 

Cr:Nd:YAG 

Shuttle  STS-72,  85 
near-Earth,  57° 

Mars  Orbiter 
Laser 

Altimeter 

Mars  laser  topography 

1064  nm 

1 

Nd:YAG 

Mars  retrograde 
-93°,  near-polar 

NEAR  Probe 

Earth-Asteroid 

rendezvous 

1064  nm 

1 

Nd:YAG 

433  Eros  proximity 
ops 

Phoenix 

Lander 

Mars  surface  probe  & 
atmospheric  LiDAR 

532,  1064  nm 

1—2  Us 
combined 

Nd:YAG 

N/A 

Shuttle 

TriDAR 

Laser  docking  & 
proximity  ops 

Unsure — 
likely  1400  or 
1550  nm 

1 

Probably 

Nd:YAG 

Shuttle  proximity 
ops  with  ISS 

ICESat/GLAS 

Clouds/aerosols, 
forest  structure, 
biomass,  ice  sheets 

532,  1064  nm 

1—2  Us 
Combined 

3  Nd:YAG 

Circular,  low- 
Earth,  near-polar 

MESSENGER 

Mercury  topography 

1064  nm 

1 

Nd:YAG 

Highly  elliptical, 
sun-shielded 

Hayabusa 

Asteroid  probe 

1064  nm 

1 

Nd:YAG 

Asteroid  rndv 

Past  Trend:  Diode-pumped,  Q-switched,  Nd:YAG  lasers;  1  -  1064  nm  beam;  planetary/lunar  orbits  less  than 
1000  km 

Present  (Systems  in  operation  at  the  time  of  writing) 

CALIPSO 

Cloud-aerosol, 

polarization 

532,  1064  nm 

1—2  Us 
Combined 

2  Nd:YAG 

Sun-synch  98°, 
low-Earth 

Lunar  Orbiter 
Laser 

Altimeter 

High  quality  lunar 
mapping  and 
topography 

1064  nm 

5 

1  Nd:YAG 
w/  split 
optics 

Elliptical,  near- 
lunar 

Future  (System  concepts  in  various  stages  of  devel 

lopment  and 

testing) 

Next  Gen. 
Integrated 
Landing 
(ALHAT) 

Spacecraft 
autonomous  landing 

Laser  Altimeter — 1570  nm  Ranges  >  2km 

Doppler  LiDAR — 1550  nm  (pre-chirp)  Ranges  <  2km 

Flash  Imaging  LiDAR  -1060  nm  Direct  hazard  avoidance 

ICESat-2/A 

TLAS 

ICES  at  follow-on 

532  nm 
(micro-pulsed) 

6 — in  3 
pairs 

Probably 

Nd:YAG 

Circular,  low- 
Earth,  near-polar 

GEDI 

Large  area  forest 
structure/biomass 

1064  nm 

14 

Nd:YAG 

ISS  attachment 

High  Altitude 
LiDAR 

Aerial  single  photon: 
wide  swath/high  res. 

532  nm 

-100 

pulselets 

Nd:YAG 

8-10  km  altitude  - 
not  yet  space  cap. 

Present/Future  Trend:  Multi-beam  (or  combination)  LiDAR  systems,  micro-pulsed  (high  Hz),  large  area 
coverage,  and  high  point  density  (spatial  resolution) 
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B. 


U.S./NASA  PLATFORMS 


1.  Beacon  Explorer — B  and  C  and  the  Goddard  Laser 

The  history  of  LiDAR  in  space  began  in  October  1964  with  ground-based  lasers 
and  satellites  carrying  cube-corner  reflectors  (Abshire  2011).  Zubritsky  explained  that, 
almost  immediately,  satellite  laser  tracking  demonstrated  a  vast  improvement  over  the 
current  method  of  tracking  with  microwave  radars;  radar  satellite  tracking  contributed 
accuracies  of  approximately  75  meters  (m)  (Zubritsky  2014).  The  same  document  makes 
an  important  comparison,  that  the  first  successful  tracking  tests  with  lasers  “reported  a 
range  accuracy  up  to  about  three  meters — about  25  times  better”  (Zubritsky  2014,  1). 

Figure  1  shows  the  Goddard  laser  (GODLAS)  with  receiving  telescope  (left)  and 
the  Beacon  Explorer-B  (BE-B)  satellite  with  its  array  of  cube-corner  reflectors  circled 
(right).  The  GODLAS  was  a  rotating-mirror,  Q- switched,  ruby  laser  mounted  on  a  NIKE- 
AJAX  radar  pedestal;  it  functioned  at  a  wavelength  of  694  nanometers  (nm)  (red)  with  20 
nanosecond  (ns)  pulses  of  approximately  0.8  joules  (J)  and  a  repetition  rate  of  1  Hertz 
(Hz)  (Degnan  2014).  Separate  elevation  and  azimuth  joysticks  moved  the  16-inch 
telescope  and  the  corresponding  955 8 A  photomultiplier  detector  (Degnan  2014). 


Figure  1.  Goddard  Laser  and  BE-B  Satellite.  Adapted  from  Abshire  (201 1). 
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The  BE-B  satellite  (Explorer  22)  was  a  52.6  kilogram  (kg)  ionospheric -research 
satellite  tasked  with  measuring  worldwide  electron  content  (NASA  Space  Science  Data 
Coordinated  Archive  [NSSDCA]  2016b).  The  NSSDCA  page  stated  that  in  addition  to 
the  laser  reflectors,  BE-B  carried  an  electrostatic  probe,  radio  beacons,  and  a  Doppler- 
navigation  system  (2016b).  This  source  confirmed  that  with  no  onboard  storage  capacity, 
data  collection  could  occur  only  when  the  satellite  was  in  view  of  a  ground  station 
(2016b).  BE-B’s  elliptical  orbit  parameters  included  an  889  kilometer  (km)  perigee, 
1,081  km  apogee,  79.7°  inclination,  and  104.8  minute  period  (NSSDCA  2016b).  Its  cube- 
comer  reflectors  returned  a  high  percentage  of  incident  light  directly  back  to  the  source 
(Degnan  2014).  BE-B’s  specifically  designed  array  of  such  reflectors  ensured  adequate 
strength  of  the  returned  signal  by  accounting  for  beam  divergence  (Degnan  2014). 

In  the  first  five  years,  the  GODLAS,  BE-B/BE-C,  other  reflector  satellites,  and 
the  deposited  lunar  retroreflectors  of  the  early  Apollo  missions  proved  the  ability  of  the 
laser  to  operate  within  the  space  environment  (Abshire  2011).  This  set  the  stage  for  a 
critical  transformation:  The  spacecraft  would  become  the  source  of  the  laser  pulse  instead 
of  the  target.  Over  subsequent  decades,  LiDAR  technology  would  expand  from  its 
humble  beginnings  into  advanced  laser  altimeters,  precision  docking  instruments,  and 
tools  for  atmospheric-constituent  measurements. 

2.  Apollo  15, 16, 17  Laser  Altimeters 

The  last  three  Apollo  missions  carried  a  flash-lamp  pumped,  mechanically  Q- 
switched,  ruby  laser  designed  by  RCA  Aerospace  (Abshire  2011).  In  the  same 
presentation,  it  is  mentioned  how  the  Apollo  laser  altimeter  came  into  existence  adapted 
from  a  tank  rangefinder  (2011).  It  weighed  23  kg,  dissipated  60  watts  of  power,  had  a 
firing  rate  of  0.05  Hz,  and  possessed  10  centimeter  (cm)  transmit  and  receive  apertures 
(Abshire  2011). 

During  Apollo  15,  16,  and  17,  the  laser  altimeter  performed  an  integral  role  in  the 
service  module’s  mission  “to  obtain  high-resolution  panoramic  and  high-quality  metric 
lunar-surface  photographs  and  altitude  data  from  lunar  orbit”  (Baldwin  1972,  1-10). 
Baldwin  announced  that  the  laser  altimeter  could  work  in  one  of  two  modes,  either  in 
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conjunction  with  or  independently  of  the  metric  camera  system  (1972).  In  the  first  mode, 
the  laser  altimeter  provided  a  corresponding  altitude  measurement  for  each  photographic 
exposure  (Baldwin  1972).  Later,  Roberson  and  Kaula  derived  an  iterative  least-squares 
triangulation  process  that  used  exposure  times,  altitude  measurements,  and  orbital  data  in 
order  to  construct  properly  georeferenced  maps  from  the  overlapping  lunar  photographs 
(Roberson  and  Kaula  1972). 

Figure  2  demonstrates  the  laser  altimeter’s  second  mode  of  operation  during 
orbits  15  and  16,  independent  of  the  cameras.  Roberson  and  Kaula  combined  altimetry 
measurements  with  orbital  data  to  create  graphs  of  the  lunar  surface  height  relative  to  a 
sphere  situated  at  the  Moon’s  center  of  mass  (1972).  The  sphere’s  radius  is  1,738  km,  and 
surface  elevation  is  represented  on  the  y-axis  in  km  above  spheroid  (Roberson  and  Kaula 
1972). 


Figure  2.  Reduced  Laser  Altimetry  for  Orbits  15/16  of  Apollo  15.  Source: 

Roberson  and  Kaula  (1972,  25-49). 


ts 

FIGURE  25-41. -Apollo  IS  b«rr  lltirnctcr  diti.(  i)  Near  tijc.  (b)  Fir  tide. 


7 


On  orbit  27,  the  Apollo  15  laser  altimeter  began  malfunctioning,  and  its 
performance  worsened  due  to  decreased  laser  output  power  until  orbit  38,  when  a 
different  anomaly  permanently  ceased  laser-altimetry  operations  (National  Aeronautics 
and  Space  Administration  [NASA]  1971).  According  to  the  Apollo  15  Mission  Report , 
NASA  never  fully  determined  the  exact  cause  of  the  degraded  output  power  anomaly,  but 
it  attributed  the  fault  to  the  laser  module  (1971).  Consequently,  future  flight  units 
included  an  automatic  power  compensation  circuit  with  a  photodiode  sensor  (1971). 
NASA  concluded  that  the  final  anomaly  likely  resulted  from  electromagnetic  interference 
due  to  a  malfunctioning  ground-te  sting  safety  relay  that  interrupted  the  proper  operation 
of  the  laser  receiver’s  video  amplifier  (1971).  As  a  result,  the  power  supply  to  the 
photomultiplier  tube  held  the  system  at  idle  voltage  for  too  long  (NASA  1971). 

Figure  3  shows  the  receiver  missing  the  return  pulse  as  a  consequence  of  the 
photomultiplier  remaining  at  idle  voltage  after  laser  firing.  NASA  would  remove  the 
safety  relay  from  future  units,  however  this  marked  the  end  of  Apollo  15  laser  altimetry 
operations  (NASA  1971). 


Figure  3.  The  Final  Apollo  15  Laser  Altimeter  Anomaly  on  Orbit  38.  Source: 

National  Aeronautics  and  Space  Administration  [NASA]  (1971,  213). 


->ooo 


o  L 


•1600  volls- 


k'ltinj  level 


7 

670  wits 


Normal  oprratit 


E 

o 

u 

c 

CL 


-2000 


lbOO  rtllS 


•1000 


-670  volts 


0  I - L 

0  2 


I  l  I _ I 

6  8  10 

T  imp ,  seconds 


Anomalous  operj 


J _ L 

12  14 


8 


The  Apollo  16  laser  altimeter  fared  somewhat  better  than  its  predecessor 
(NSSDCA  2016a).  Per  NSSDCA,  the  system’s  average  reliability  had  decreased  by  orbit 
60  to  approximately  50%  (2016a).  Unfortunately,  the  laser  altimeter  failed  three  orbits 
later  during  the  final  run  of  photograph  collection  (NSSDCA  2016a). 

The  Apollo  17  laser  altimeter  operated  successfully  for  the  duration  of  the  mission 
(Kaula  et  al.  1974).  Kaula  et  al.  stated  that  concern  about  the  prior  failures  caused  the 
astronauts  of  Apollo  17  to  avoid  excessive  use  of  the  laser  altimeter  early  in  the  mission 
to  ensure  its  later  availability  to  support  photographic-mapping  operations  (1974).  Figure 
4  indicates  the  position  of  the  laser  altimeter  on  the  Apollo  17  Scientific  Instrument 
Module.  The  photograph  is  from  the  Command  Module  in  lunar  orbit. 


Figure  4.  Apollo  17  Scientific  Instrument  Module  and  Layout  Diagram. 

Adapted  from  Lunar  and  Planetary  Institute  and  Universities  Space 
Research  Association  (2016). 
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Kaula  et  al.  combined  laser  altimetry  data  from  the  Apollo  15,  16,  and  17 

missions  to  derive  average  elevations  for  important  lunar  terrains  in  relation  to  the  1,738 

km  (mean  radius)  spheroid:  “farside  terrae  +1.8  km,  nearside  terrae  -1.4  km,  ringed  maria 

-4.0  km,  and  other  maria  -2.3  km”  (Kaula  et  al.  1974,  3049).  Their  analyses  of  these 

laser-derived  elevation  differences,  along  with  gravitational  and  surface  density 
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measurements,  led  them  to  conclusions  about  the  composition  and  density  of  the  lunar 
core:  “It  is  unlikely  that  any  iron  core  is  greater  than  400  km  in  radius.  A  more  likely 
model  is  a  gradual  increase  in  density  from  3.33  g/cm-3  to  3.46  g/cm-3  at  1,000  km 
depth,  compatible  with  Mg/Fe  ratios  plausible  for  mare  basalt  parent  material”  (Kaula  et 
al.  1974,  3049). 

3.  Clementine 

Launched  25  January  1994,  the  Clementine  lunar  probe  originated  from  a  NASA/ 
Ballistic  Missile  Defense  Organization  partnership.  According  to  the  National  Space 
Science  Data  Center  (NSSDC),  Clementine’s  mission  was  “to  test  sensors  and  spacecraft 
components  under  extended  exposure... and  to  observe  the  Moon  and  near- Earth  asteroid 
1620  Geographos”  (National  Space  Science  Data  Center  [NSSDC]  1994,  1).  NSSDC 
claimed  that  after  completing  two  months  of  lunar  mapping  operations  and  an  Earth 
transfer  orbit,  Clementine  suffered  from  an  attitude  control  thruster  malfunction  that 
exhausted  the  thruster’s  remaining  hydrazine  fuel  supply  and  spun  the  spacecraft  at  80 
revolutions  per  minute  (1994).  Therefore,  NASA  scrapped  the  asteroid  rendezvous 
portion  of  the  mission,  and  Clementine  headed  into  the  Van  Allen  radiation  belts  to  test 
its  remaining  functional  components  (1994).  The  mission  ended  in  June  1994  when 
intelligible  transmissions  from  Clementine  to  Earth  ceased  (NSSDC  1994). 

Figure  5  provides  a  drawing  of  the  1.88  m  by  1.14  m  (main-body  dimensions), 
227  kg  Clementine  spacecraft  and  its  diode-pumped  chromium  neodymium-doped 
yttrium  aluminum  garnet  (Cr:Nd:YAG)  solid-state  compact  LiDAR  (NSSDC  1994).  The 
NSSDC  details  page  informed  us  that  the  Clementine  LiDAR  functioned  at  532  nm  and 
1064  nm  wavelengths  (green  and  infrared)  with  pulse  energies  of  nine  and  180  milli- 
joules  (mJ),  respectively  (1994).  Both  the  infrared  ranging  and  green  active-imaging 
wavelengths  had  a  pulse-width  of  less  than  10  ns  (NSSDC  1994).  During  lunar-mapping, 
Clementine’s  LiDAR  frequently  outperformed  its  maximum  effective  range  of  500  km 
(NSSDC  1994).  The  Clementine  project  represented  a  vast  improvement  over  the  Apollo 
laser  altimeters.  Researchers  geolocated  72,300  of  Clementine’s  600,000  total  laser  shots, 
and  accuracy  averaged  90  m  vertically  and  3  km  horizontally  (Neumann  2001). 
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Figure  5.  Drawing  of  the  Clementine  Spacecraft  and  its  LiDAR  System. 

Source:  Williams  (2011). 


Figure  6  is  an  equal-area  projection  topographic  map  of  the  Moon  derived  from 
Clementine  LiDAR  data.  Elevations  range  from  about  -8  km  (purple)  to  +4  km  (pink/ 
white).  Also,  the  Lunar  and  Planetary  Institute  informs  us  that  the  farside  image  is  tilted 
in  order  to  show  the  low-lying,  2,500  km  wide,  Aitken  Basin  near  the  lunar  south  pole 
(Lunar  and  Planetary  Institute  2016). 


Ligure  6.  Clementine  Topographic  Map  of  the  Moon.  Source:  Lunar  and 

Planetary  Institute  (2016). 
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4.  Space  Shuttle  LiDAR  In-space  Technology  Experiment  and  Laser 
Altimeter  Experiments 

In  1994,  NASA  released  a  factsheet  that  contained  an  interesting  eye-safety  laser 
disclaimer  and  a  warning  to  the  international  science  community  to  avoid  using  a 
medium-plus  sized  telescope  to  view  the  Space  Shuttle  Discovery  during  STS-64 
(McCormick  et  al.  1994).  The  Space  Shuttle  LiDAR  In-space  Technology  Experiment 
(LITE)  was  a  first:  an  atmospheric  LiDAR  survey  conducted  by  an  Earth-orbiting 
spacecraft  (Abshire  2011).  Figure  7  provides  an  illustrated  diagram  of  the  payload  bay 
apparatus  and  the  LITE  mission  concept. 

According  to  Abshire,  Discovery  maintained  a  circular  low-Earth  orbit  at  206  km 
for  nine  days  of  data  collection  (2011).  He  explained  how  the  thin  laser  pulses  traveled 
downward  through  the  mid- stratosphere  and  the  troposphere,  diverging  into  300  m 
ground  spots  spaced  740  m  along-track  (2011).  Overall,  LITE  collected  53  hours  of  near¬ 
nadir  atmospheric  backscatter  profiles  on  the  order  of  10  seconds  each  (Abshire  201 1). 


Figure  7.  LITE  Sample  Mission  and  Laser  Apparatus  Cartoons.  Adapted  from 

McCormick  et  al.  (1994). 
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Altogether,  the  LITE  transmitter  apparatus,  boresight  assembly,  telescope 
receiver,  camera,  and  system  electronics  weighted  1,304  kg  and  consumed  3,030  watts 
(McCormick  et  al.  1994).  McCormick  commented  that  the  pressurized  and  thermal- 
controlled  laser  transmitter  module  encased  two  identical,  flash-lamp  pumped,  Q- 
switched  Nd:YAG  lasers  (1994).  He  added  that  the  two  lasers  never  fired  simultaneously 
to  ensure  redundancy,  and  the  initial  laser  pulse  created  was  a  1064  nm  (infrared)  beam 
(1994).  McCormick  went  on  to  say  that  this  pulse  was  temporarily  split  and  then 
recombined  into  one  beam  containing  three  wavelengths — harmonic  generators  converted 
two  of  the  three  split-beams  into  the  additional  532  nm  (green)  and  355  nm  (ultraviolet) 
wavelengths  (1994).  Additionally,  LITE  used  gimballed  turning  prisms  to  point  the  laser 
boresight  at  the  Earth,  and  the  receiver  encased  a  Ritchey-Chretien  Cassegrain  telescope 
with  beryllium  as  the  primary  mirror  and  quartz  as  the  secondary  mirror  (McCormick  et 
al.  1994).  The  aft  optics  assembly  contained  photomultiplier  tubes  for  detecting  the 
ultraviolet  and  green  ground/atmospheric  returns  and  a  silicon  avalanche  photodiode  for 
the  infrared  signal  (McCormick  et  al.  1994). 

Almost  10  years  after  program  commencement,  LITE  accomplished  the  goals  of 
its  creators  and  the  LITE  Science  Steering  Group  “to  demonstrate  operation  of  a  LiDAR 
in  space  and  the  maturation  of  LiDAR  technology  toward  the  development  of  future 
systems  on  free-flying  satellite  platforms”  (Winker,  Couch,  and  McCormick  1996,  1). 
Winker,  Couch,  and  McCormick  claimed  that  LITE  would  pave  the  way  for  future 
longer-term  atmospheric  and  climate-monitoring  satellites,  and  they  were  right:  ICESat, 
CALIPSO,  and  the  upcoming  ICESat-2  (Winker,  Couch,  and  McCormick  1996). 

Within  three  years  of  LITE,  shuttle-based  LiDAR  would  conduct  two  additional 
missions,  albeit  for  a  traditional  topographic  purpose.  Shuttle  Laser  Altimetry  (SLA) 
Missions  1  and  2  “demonstrated  the  effectiveness  of  orbital  laser  altimeters  for  terrestrial 
geodesy  despite  the  inconvenience  of  using  the  shuttles  as  an  orbital  platform”  (Neumann 
2001,  73).  Lurthermore,  in  the  publication  “Some  Aspects  of  Processing  Extraterrestrial 
LiDAR  Data — Clementine,  NEAR,  MOLA,”  Neumann  mentioned  that  left-over  spares 
from  the  Mars  Orbiter  Laser  Altimeter  (MOLA)  project  plus  additional  components  gave 
rise  to  the  SLA  (2001).  Although,  Neumann  did  not  disclose  much  further,  he  did  provide 
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some  statistics  (listed  in  Table  2)  on  the  MOLA-offshoot  that  became  the  LiDAR  system 
for  the  SLA-01  and  SLA-02  missions. 


Table  2.  SLA  Mission  Quick  Stats.  Adapted  from  Neumann  (2001). 


Mission 

STS 

Date 

Laser 

Shots 

Fired 

Shots 

Geolocated 

Horizontal 

Accuracy 

Vertical 

Accuracy 

SLA-01 

STS-72 

Jan  1996 

Cr:Nd:YAG 
at  10  Hz 

3  Million 

1,203,000 

40  m 

2.78  m 

SLA-02 

STS-85 

Aug  1997 

Cr:Nd:YAG 
at  10  Hz 

3  Million 

2,090,000 

40  m 

6.74  m 

Unlike  the  atmospheric-sensing  LITE,  SLA  only  needed  one  40  mJ,  1064  nm 
(infrared)  laser  beam  (Bufton  et  al.  1995).  Bufton  et  al.  described  that  the  SLA  pulse- 
width  was  10  ns,  and  the  resulting  terrestrial  footprint  was  100  m  in  diameter  (1995).  A 
silicon  avalanche  diode  detected  the  echoed  waveforms  (1995).  Next,  they  stated  that  a 
pulse  waveform  digitalizer  interpreted  the  return  signal  in  order  to  resolve  the  terrain 
topography  and  vegetation  height-above-ground  profiles  within  each  laser  footprint 
(1995).  Furthermore,  reducing  the  waveform  resolution  and  sampling  for  a  longer 
duration  enabled  the  SLA  to  function  as  an  atmospheric  and  cloud-aerosol  LiDAR 
(Bufton  et  al.  1995).  Attached  as  a  Hitchhiker  payload,  the  SLA  transmitted  its  pulses 
directly  along  nadir  while  the  Shuttle  orbited  at  around  300  km;  along-track  sampling 
was  similar  to  LITE,  approximately  740-750  m  (Sun  et  al.  2003). 

Figure  8  shows  a  composite  map  of  the  mostly  tropical  (28.45°  N  and  S)  SLA-01 
ground  tracks  and  the  57°  inclination  SLA-02  tracks.  Harding  et  al.  used  data  collected 
from  the  SLA-01  mission  over  “Africa,  southern  Asia,  central  South  America,  and 
Australia”  to  test  the  vertical  accuracy  of  the  U.S.  Geological  Survey’s  Global  30  Arc- 
second  Digital  Elevation  Model  (DEM)  Dataset  (GTOPO30),  a  DEM  with  30  arc-second 
(~1  km)  grid  spacing  (Harding  et  al.  1998,  1).  For  example,  Australia  had  29,139 
discrepancies  between  SLA-01  measurements  and  the  GTOPO30,  with  an  average  error 
of  -21.72  m  and  a  standard  deviation  of  48.92  m;  Africa’s  discrepancies  totaled  244,640, 
with  a  mean  error  of  -1.40  m  and  a  standard  deviation  of  44.75  m  (1998).  Additionally, 
Harding  et  al.  realized  that  the  SLA’s  ability  to  achieve  an  accurate  ground  return  in  a 
forested  area  depended  significantly  on  the  density  of  the  vegetation  (1998).  Also,  they 
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noted  that  the  SLA  mapped  only  the  cloud-top  layer  when  it  operated  above  opaque 
clouds  (Harding  et  al.  1998). 


Figure  8.  Earth  Surface  Coverage  for  Both  Shuttle  Laser  Altimetry  Missions. 

Source:  Carabajal  et  al.  (2000,  2). 
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Behn  and  Zuber  combined  SLA-01  LiDAR  data  and  TOPEX/POSEIDON  radar 
altimetry  to  study  ocean  topography  (Behn  and  Zuber  2000).  They  discovered  a  2.33  m 
difference  (root  mean  square)  between  the  University  of  Texas  Mean  Sea  Surface  Model 
(CSRMSS95)  and  the  more  than  one  million  SLA-01  sea  surface  measurements  (2000). 
Moreover,  Figure  9  illustrates  that  the  SLA  LiDAR  could  successfully  maintain  track 
even  when  crossing  the  land-sea  threshold  in  coastal  areas  with  drastically  varying 
elevations  (Behn  and  Zuber  2000).  The  top  region  presented  in  Figure  9  is  a  selection  of 
the  South  American  coast  and  the  bottom  region  is  a  segment  of  the  Red  Sea.  Also,  the 
laser  points  are  represented  by  red  dots,  and  the  land  outline  is  a  conjunction  of  National 
Geophysical  Data  Center  Terrain  Base  elevations  and  over-ocean  segments  of  the  Earth 
Gravitational  Model  1996  (EGM96)  geoid  (Behn  and  Zuber  2000). 

Sun  et  al.  compared  SLA-02  altimetry  data  to  the  corresponding  surface  heights 
from  the  year  2000  Shuttle  Radar  Topography  Mission  (SRTM)  to  create  a  highly 
accurate  and  digital  global-topographic  database  (Sun  et  al.  2003).  Their  first  step  was 
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filtering  out  cloud  returns  by  eliminating  any  point  with  an  elevation  above  3  km  (2003). 
Their  follow-on  tasks  involved  analyzing  all  waveforms  and  terrain  returns,  then 
associating  them  to  Landsat  imagery  (2003).  The  LiDAR  data  was  instrumental  in 
validating  the  radar  data  and  led  to  their  strong  conclusion:  “The  SRTM  DEM  proves  to 
be  the  best-known  DEM  ever  generated  at  the  global  scale  with  its  consistency  and 
overall  accuracy”  (Sun  et  al.  2003,  410). 


Figure  9.  Two  SLA-01  Profiles  Illustrating  the  Ability  to  Track  Across  the 
Land-Sea  Interface.  Source:  Behn  and  Zuber  (2000,  1436). 


Only  two  of  the  four  originally  planned  SLA  missions  occurred,  nevertheless 

Bufton  et  al.  correctly  predicted  that  the  SLA  would  serve  as  an  important  precursor  for 

future  space-based  LiDAR  systems,  including  the  Geoscience  Laser  Altimeter  System 

(GLAS) — eventually  to  fly  onboard  ICESat  (Bufton  et  al.  1995).  Indeed,  LITE  and  SLA 
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marked  the  beginning  of  an  ongoing  NASA  Goddard  endeavor  to  transcend  aircraft 
LiDAR  systems  into  “space-based  laser  remote  sensing  devices”  (Bufton  et  al.  1995,  1). 

5.  Mars  Global  Surveyor/Mars  Orbiter  Laser  Altimeter 

Over  the  course  of  its  decade-long  mission  from  November  1996-2006,  the  Mars 
Global  Surveyor  (MGS)  spacecraft  and  its  LiDAR,  the  Mars  Orbiter  Laser  Altimeter 
(MOLA),  provided  laser  topographic  mapping  across  the  entire  surface  of  the  red  planet 
(NASA  2007b).  Figure  10  offers  a  rendition  of  the  MGS  and  a  photograph  of  MOLA. 


Figure  10.  Mars  Global  Surveyor  and  MOLA.  Left  Source:  NASA  Jet  Propulsion 
Laboratory  [JPL]/California  Institute  of  Technology  [Caltech]  (2012). 
Right  Source:  Smith,  David,  and  Zuber  (2007). 


According  to  Smith,  the  25.9  kg  and  30.9  watt  MOLA  was  nearly  identical  to  the 
Mars  Observer  Laser  Altimeter  (same  acronym)  that  had  flown  onboard  the  ill-fated  Mars 
Observer  Probe  three  years  prior  (Smith  1996).  He  mentioned  that  MOLA’s  non- 
topographic  objectives  included  measuring  the  vertical  profiles  of  carbon  dioxide  and 
water  clouds  in  the  Martian  atmosphere  (1996).  Furthermore,  Smith  and  other  scientists 
combined  MOLA  LiDAR  data  with  measurements  from  the  spacecraft’s  other 
instruments,  such  as  the  thermal  emission  spectrometer  and  the  magnetometer,  to  study 
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the  internal  compositions  of  the  two  polar  caps  and  the  other  planetary  structures  on  Mars 
(Smith  1996). 

Figure  11  shows  a  diagram  of  the  laser  waveform  produced  by  MOLA’s  Q- 
switched,  diode-pumped,  Nd:YAG  laser  that  operated  at  10  Hz,  40  to  45  mJ,  8  ns  pulse- 
widths,  and  a  1064  nm  (infrared)  wavelength  (Abshire  2011).  Abshire  stated  that  the 
silicon  avalanche  photodiode  analog  detector  was  effective  at  catching  a  minimum  return 
signal  of  230  photons  (represented  by  the  line  labelled  detection  threshold);  pulse  round- 
trip  travel  time  and  the  corresponding  laser-derived  ranges  enabled  the  calculation  of 
surface  elevation  (2011).  The  shape  of  the  returned  laser  waveform  described  terrain 
roughness,  and  a  comparison  of  the  transmitted  versus  received  pulse  energies  yielded  a 
means  to  predict  surface  reflectance  (Abshire  2011).  MOLA’s  horizontal  and  absolute 
vertical  resolutions  were  160  m  and  30  m,  respectively  (Smith  1996). 


Figure  11.  Sample  MOLA  Measurement.  Source:  Abshire  (2011,  11). 
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An  example  of  one  of  about  67 1  million  MOLA  laser  shots. 


The  MGS  mission  represents  a  historic  success.  MOLA  data  vastly  improved  the 
scientific  community’s  overall  understanding  of  Mars  (NASA  JPL/Caltech  2016).  Prior 
knowledge  about  the  red  planet  came  almost  exclusively  from  less  accurate  and  outdated 
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Mariner,  Viking,  and  Earth-based  radar  measurements  (Smith  et  al.  2001).  Per  Smith  et 
al.,  after  spending  the  first  two  years  on  station  in  capture  and  aerobrake  testing  orbits, 
the  MGS  spacecraft  commenced  its  LiDAR-mapping  mission  in  February  of  1999  from  a 
400  km  circular,  near-polar,  orbit  at  92.8°  inclination  (2001).  Figure  12  provides 
topographic  maps  of  the  Martian  north  and  south  poles  (elevations  in  m).  Figure  13  is  a 
side-by-side  topographic  map  of  the  eastern  and  western  hemispheres  (elevations  in  km). 


Figure  12.  MOFA  FiDAR  Elevation  Maps  of  Mars’  North  and  South  Poles. 

Source:  NASA  JPF/Caltech  (2016). 

180'  O' 


180' 


BOOO  -4000  0  4000 

Topography  [m] 

North  pole  (Left)  labels:  Olympia  Planum — Icecap  plateau  at  (180°  E,  80°  N).  Boreale, 
Chasma — The  deep  trench  in  the  ice  cap  at  (330°  E,  85°  N). 


BOOO 
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southern  icecap.  Australe,  Chasma — A  canyon  near  (90°  E,  85°  S). 
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Figure  13.  MOLA  LiDAR  Elevation  Map  of  Mars’  Eastern  and  Western 
Hemispheres.  Source:  NASA  JPL/Caltech  (2016). 


-8 

•4 

0  4 

8 

12 

Altitude  [km] 

6.  Near  Earth  Asteroid  Rendezvous  Spacecraft — Asteroid  433  Eros 

From  2000  to  2001,  NASA’s  Near  Earth  Asteroid  Rendezvous  (NEAR) 
spacecraft,  nicknamed  Shoemaker,  performed  a  scientific  study  of  the  asteroid  433  Eros 
from  a  predominantly  35  km  orbit  (Cole  1998).  Cole  explained  that  the  NEAR  Laser 
Rangefinder’s  (NLR)  mission  was  to  collect  precise  laser  altimetry  measurements  that 
would  supply  scientists  with  high  quality  data  products  unobtainable  from  Earth-based 
asteroid  research  (1998).  Also,  he  revealed  that  in  order  to  carry  out  its  extended  deep- 
space  tasking,  NEAR  maintained  orbital  trajectories  which  kept  its  LiDAR  and  other 
instruments  directed  at  433  Eros  while  simultaneously  pointing  the  topside 
communications  antenna  and  solar  arrays  at  the  Earth  and  Sun,  respectively  (1998). 
Figure  14  demonstrates  this  feat  of  spacecraft  attitude  control,  and  it  also  shows  a  single 
NLR  footprint  of  12  m — not  to  scale  on  the  estimated  36  km  x  15  km  x  13  km  (length, 
width,  height)  asteroid. 
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Figure  14.  Drawing  of  Orbital  Requirements  for  NEAR  and  Photograph  of  Its 
Laser  Rangefinder.  Source:  Cole  (1998,  143,  145). 


Figure  14  also  provides  a  picture  of  the  LiDAR  device  that  successfully  mapped 
433  Eros  from  ranges  up  to  100  km  (Cole  1998).  Cole  commented  that  pre-launch  testing 
of  the  NLR  yielded  an  average  range  accuracy  of  32  cm  and  a  beam  divergence  of  235 
micro-radians  (Cole  1998).  Per  Abshire,  the  4.9  kg  NLR  operated  at  1  Hz,  16.5  watts  of 
power  (average),  15  mJ  pulse  energies,  12  to  15  ns  pulse-widths,  and  a  wavelength  of 
1064  nm  (infrared)  (Abshire  2011).  Also,  the  laser  transmitter  was  a  Q-switched,  diode- 
pumped,  Nd:YAG  system  (Abshire  2011).  Cole  mentioned  that  the  NLR’s  avalanche 
photodiode  detector  employed  advanced  timing  electronics  for  time-of-flight  (TOF) 
calibrations  and  for  accurate  detection  of  the  laser  backscatter  from  the  asteroid’s  surface 
(Cole  1998). 

In  their  report  titled  “The  Shape  of  433  Eros  from  the  NEAR- Shoemaker  Laser 
Rangefinder,”  Zuber  et  al.  presented  topographic  findings  and  diagrams  based  on  the 
NLR-collected  data  (Zuber  et  al.  2000).  Figure  15  is  a  radius  map  of  the  asteroid  using  a 
cylindrical  projection;  longitude  is  positive  in  the  easterly  direction  from  0°  to  360°,  and 
latitude  ranges  from  -90°  to  90°.  Zuber  et  al.  indicated  that  latitude  separation  is 
approximately  54  m  at  the  asteroid’s  poles,  expanding  to  307  m  at  the  equator  (2000).  Of 
note,  the  low  regions  named  Himeros  and  Psyche  center  themselves  at  75°  East  and  270° 
East,  respectively. 
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Figure  15.  Radius  Map  of  433  Eros  in  Cylindrical  Projection. 

Source:  Zuber  et  al.  (2000,  2098). 
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Furthermore,  the  NLR  enabled  Zuber  et  al.  to  determine  the  overall  shape  and 
slopes  of  433  Eros  (2000).  They  concluded  that  the  asteroid  generally  slopes  between  five 
and  15  degrees  with  exception  of  the  two  sharply  concave  regions  Himeros  and  Psyche 
that  together  form  a  wide  depression  running  north  to  south,  slightly  offset  toward  the 
positive  x-direction  (Zuber  et  al.  2000).  Figure  16  provides  two  3D  views  of  433  Eros 
with  slopes  color-scaled  from  a  0°  low  (purple)  to  a  40°  high  (white).  The  left  view  is 
from  (300°E,  30°S),  and  the  right  view  is  from  (0°E,  30°S). 
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Figure  16.  Two  Views  of  a  3D  Shape  Model  of  433  Eros  from  the  NLR. 

Source:  Zuber  et  al.  (2000,  2099). 
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7.  Phoenix  Mars  Lander 

Landing  on  the  red  planet  in  May  2008,  the  mission  of  Phoenix  Mars  Lander  was 
to  investigate  a  region  on  Mars  where  scientists  expected  the  topsoil  to  contain  a 
significant  amount  of  frozen  water  (NASA  2008). 1  NASA  explained  that  Phoenix  would 
conduct  ground-based  weather  research  near  the  Martian  north  pole,  specifically  focusing 
on  the  seasonal  water  cycle  where  ground-ice  sublimes  into  atmospheric  water  vapor  and 
later  depositions  back  into  solid  form  (2008).  Phoenix’s  meteorological  LiDAR  measured 
the  downward-reflected  laser  backscatter  in  order  to  collect  altitude,  abundancy,  and 
sizing  data  for  particles  of  dust  and  ice  in  the  atmosphere  (NASA  2008). 

Ligure  17  presents  an  artistic  rendition  of  the  probe’s  robotic  arm  hard  at  work 
collecting  soil  samples  for  onboard  testing.  Meanwhile,  its  top-mounted  LiDAR  conducts 
laser  measurements  of  the  Martian  atmosphere.  Ligure  17  also  displays  a  true-color 
surface  stereo  photograph  taken  by  Phoenix  on  25  June  2008. 


1  Earlier  gamma-ray  spectrometer  observations  by  the  Mars  Odyssey  Orbiter  aided  in  selection  of  the 
Phoenix  landing  site;  Odyssey  measurements  suggested  that  permafrost  soil  existed  near  the  planet’s  poles 
(NASA,  2008). 
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Figure  17.  The  Phoenix  Mars  Lander.  Source:  NASA  (2015). 


The  Canadian  Space  Agency  contributed  the  Phoenix  weather  instrument  suite, 
and  the  LiDAR  came  courtesy  of  MDA  Space  Systems,  Ontario  and  Optech  Inc., 
Toronto.  Whiteway  et  al.  remarked  that  the  laser  transmitted  at  100  Hz  via  a  passively  Q- 
s witched,  diode -pumped,  Nd:YAG  assembly  that  produced  a  pulse- width  of  10  ns  and  a 
total  pulse  energy  of  0.7  mJ:  0.3  mJ  at  1064  nm  (infrared)  and  0.4  mJ  at  532  nm  (green) 
(Whiteway  et  al.  2008).  They  stated  that  a  potassium- titanyl  phosphate  crystal  created  the 
visible  green  wavelength  from  part  of  the  original  infrared  beam  through  frequency 
doubling  (2008).  Additionally,  the  10  cm  diameter  receiver  employed  a  silicon  avalanche 
photodiode  detector  that  recorded  the  return  signals  via  14-bit  analog  and  photon¬ 
counting  methods  (Whiteway  et  al.  2008). 

Overall,  Phoenix  exceeded  expectations  by  significantly  outlasting  its  intended 
mission  duration  of  90  Martian  solar  days  (sols)  (Whiteway  et  al.  2008).  A  sol  is 
approximately  24  hr  39  min  35  sec.  Figure  18  provides  a  LiDAR-derived  vertical  profile, 
local  time  versus  altitude,  of  the  Martian  clouds  directly  above  the  probe  on  sol  99  (3 
September  2008);  starting  at  0500  local,  ice  crystals  are  detected  falling  from  the  clouds 
at  an  altitude  of  about  3.25  km  (NASA/JPL-Caltech  2016). 


24 


Figure  18.  Phoenix  LiDAR  Profile  of  Martian  Clouds  with  Falling  Ice  Crystals. 

Source:  NASA/JPL-Caltech  (2016). 
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Collected  on  3  September  2008. 


8.  Space  Shuttle  Triangulation  +  LiDAR  Automated  Rendezvous  and 
Docking 

Although  the  purpose  of  the  LiDAR  system  discussed  here  differs  substantially 
from  nearly  all  of  the  others,  it  is  nonetheless  worthy  of  inclusion  as  a  two-fold  testament 
to  the  future  of  precision  spacecraft  rendezvous  and  to  the  honored  memory  of  the  Space 
Shuttle.  Figure  19  is  a  photograph  taken  from  the  International  Space  Station  (ISS)  of  the 
Space  Shuttle  Discovery  testing  the  Triangulation  +  LiDAR  (TriDAR)  spacecraft 
precision  docking  system  during  STS- 131  in  April  2010  (NASA  Human  Space  Flight 
2014).  Neptec  Design  Groups’  TriDAR  is  a  self-sufficient  “vision”  laser  system  that 
applies  its  own  models  and  protocols  during  proximity  operations  without  the  assistance 
of  retroreflectors  or  other  inputs  from  the  rendezvous  target  (Ruel,  Luu,  and  Berube  2010, 
15).  Representing  another  international  endeavor  between  NASA  and  the  Canadian  Space 
Agency,  the  TriDAR  completed  three  space-testing  flights  on  STS-128,  STS-131,  and 
STS-135  (Christian  and  Cryan  2013). 
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Figure  19.  Discovery  Docking  with  the  ISS  Using  TriDAR.  Adapted  from  NASA 

Fluman  Space  Flight  (2014). 


STS-131,  April  2010.  TriDAR  is  circled. 


On  these  STS  missions,  TriDAR  demonstrated  the  ability  to  derive  the  target 
spacecraft’s  shape,  proximity,  and  relative  velocity  from  its  analyses  of  successive  3D 
point  clouds  (Christian  and  Cryan  2013).  Furthermore,  starting  with  STS-131,  TriDAR 
successfully  maintained  a  real-time  track  on  the  ISS  during  a  360-degree 
circumnavigation  without  any  external  automation  or  human-controlled  inputs  (Ruel, 
Luu,  and  Berube  2010). 

Figure  20  graphs  the  three-axis  docking  approach  of  Discovery  with  the  ISS 
during  the  first  TriDAR  space-test  in  August  2009.  Relative  distances  on  all  y-axes  are  in 
feet,  and  each  graph  includes  three  vertical  lines  to  distinguish  among  the  final  approach, 
triangulation  transition,  and  capture  phases.  Ruel,  Luu,  and  Berube  noted  that  the  total 
process  took  approximately  30  minutes  from  the  go-for-dock  command  (2010). 
Moreover,  the  TriDAR  system  continuously  fed  the  Space  Shuttle’s  pilots  highly  accurate 
closing  data,  including  automatically  updated  state  vectors  giving  the  relative  position  of 
the  ISS  throughout  docking  and  undocking  operations  (Ruel,  Luu,  and  Berube  2010). 
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Figure  20.  Three-Axis  Docking  Distances  from  Discovery  (STS-128)  to  ISS. 

Source:  Ruel,  Luu,  and  Berube  (2010,  18). 
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According  to  Neptec,  the  TriDAR  system  uses  laser-triangulation  methods  similar 
to  those  applied  by  the  Space  Shuttle’s  Laser  Camera  System  (Ruel,  Luu,  and  Berube 
2010).  Ruel,  Luu,  and  Berube  mentioned  that  a  single  laser  optic  conducted  both 
triangulation  and  LiDAR  functions  to  collectively  improve  the  system’s  capability  to 
accurately  detect  objects  in  space  (2010).  Additionally,  the  system  employs  a  thermal 
detector  to  extend  initial  target  detection  and  subsequent  rendezvous  guidance  beyond  the 
maximum  range  of  the  laser  component  alone  (a  few  km)  (2010).  With  the  thermal 
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imager  operational,  the  TriDAR  system  successfully  detected  the  ISS  at  a  range  of  just 
over  43  km  (Ruel,  Luu,  and  Berube  2010). 

Lastly,  Figure  21  is  a  3D  TriDAR- generated  laser  image  of  the  ISS.  Space  Shuttle 
Discovery  compiled  the  image  during  undocking  operations  and  throughout  the  shuttle’s 
subsequent  circumnavigation  of  the  station  (Ruel,  Luu,  and  Berube  2010).  The  image  is 
colored  by  relative  ranges  from  near  to  far:  red,  yellow,  green,  and  blue.  Despite 
Discovery  passing  behind  the  ISS  while  mapping,  the  opposite  side  is  still  represented  in 
the  TriDAR  point  cloud  as  farther-away  (blue)  due  to  the  outward  spiral-type  maneuver 
taken  by  the  Shuttle  around  the  station  (Ruel,  Luu,  and  Berube  2010). 


Figure  21.  TriDAR  3D  Mapping  of  ISS  During  STS-128  Undocking  Operations. 

Source:  Ruel,  Luu,  and  Berube  (2010,  19). 


9.  Ice,  Cloud,  and  Land  Elevation  Satellite/Geoscience  Laser  Altimeter 
System 

In  the  early  2000s,  the  Shuttle’s  LITE  and  Laser  Altimetry  missions  of  the  1990s 
became  surpassed  in  duration,  scope,  and  magnitude.  NASA’s  Ice,  Cloud,  and  Land 
Elevation  Satellite  (ICES at)  carried  one  scientific  instrument,  the  Geoscience  Laser 
Altimeter  System  (GLAS);  the  final  stats  are  seven  years  in  operation  (2003-2010)  and 
1.98  billion  laser  shots  (Abshire  2011).  ICESat  had  numerous  objectives,  but  the 
principle  ones  involved  measuring  the  changing  elevations  and  mass  balances  of  the  ice 
sheets  and  landmasses  located  at  the  Earth’s  Arctic  and  Antarctic  regions  (NSSDC  2003). 
NASA  scientists  would  use  this  data  to  calculate  the  current  and  estimated-future  effects 
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of  climate  change  on  global  sea  levels  and  polar  sea-ice  volumes  (2003).  ICESat’s 
secondary  missions  included  additional  LiDAR- specialized  tasks,  such  as  measuring  the 
vertical  structures  and  heights  of  clouds  and  aerosols;  mapping  terrain  topography, 
reflectivity,  and  roughness  (outside  of  the  polar  regions);  and  analyzing  forest  structures 
and  biomass,  vegetation  heights,  and  snow  patterns  (NSSDC  2003).  Figure  22  provides  a 
rendition  of  ICESat  and  a  photograph  of  GLAS  under  construction. 


Figure  22.  Artist  Rendition  of  ICESat  on  Orbit  and  Boresight  View  of  GLAS. 

Adapted  from  NASA  Goddard  (2015). 


After  initial  on-orbit  testing,  ICESat  repositioned  into  its  94°  low-earth,  nearly 
circular,  orbit  with  an  apogee  of  594  km  and  a  perigee  of  586  km  (NSSDC  2003).  Figure 
23  demonstrates  ICESat  using  Global  Positioning  System  (GPS)  receivers  to  maintain 
accurate  orbital  timing  and  positioning.  Additionally,  it  shows  ICESat  employing  star 
trackers  for  precision  spacecraft  attitude  control  and  laser  boresight  pointing.  The  ground 
tracks,  transmitted  along  nadir,  are  spaced  170  m  apart  with  each  laser  footprint 
measuring  70  m  in  diameter  (Abshire  201 1). 
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Figure  23.  GLAS  Measurement  Approach.  Source:  Abshire  (2011,  24). 


To  achieve  the  intended  three-  to  five-  year  data  collection  period,  GLAS  carried 
three  transmitter  modules,  and  each  one  had  a  40  Hz  Nd:YAG  laser  that  produced  a 
combined  1064  nm  (infrared)  and  532  nm  (green)  pulse  (Abshire  2011).  Abshire  stated 
that  pulse-widths  were  4  ns  in  duration  (2011).  Additionally,  he  specified  that  a  silicon 
avalanche  photodiode  served  as  the  infrared  detector  capable  of  recording  the  return 
signal  at  1  Gigahertz  (GHz)  (2011).  Photon-counting  detectors  captured  the  green 
ground-return  pulses  and  atmospheric  backscatter  (Abshire  2011).  The  infrared 
component  measured  surface  altimetry  and  thick-cloud  elevation  data  at  accuracies  of  15 
cm,  and  the  green  component  measured  thin-cloud  and  atmospheric-aerosol  distributions 
(Zwally  et  al.  2002). 

The  literature -base  of  research  conducted  from  ICESat/GLAS  measurements  is  as 
extensive  as  the  numerous  raw  and  processed  data  products  themselves.  The  author 
discovered  more  scientific  publications  hereunto  pertaining  than  for  any  of  the  other 
systems  discussed  in  this  review.  Hieu  van  Duong  analyzed  the  quality  of  large-footprint 
full  waveforms  and  elevation  data  for  selected  areas  in  Europe,  discovering  that  the 
available  ICESat  data  conformed  reasonably  well  to  a  Gaussian  decomposition  approach 
for  waveform  analysis  despite  the  low  number  of  identical-repeat  footprints  across  the 
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European  continent  (Duong  2010).  Likewise,  Yadav’s  study  of  the  complex,  multi- 
peaked,  ICES  at  waveforms  from  ravines  and  other  areas  of  rough  terrain  created 
improved  DEMs  for  those  areas  as  compared  to  DEMs  generated  exclusively  from  the 
Earth  Observing  Cartosat  satellites  (Yadav  2010). 

Figure  24  provides  examples  of  Duong’s  Gaussian  decomposition  approach  for 
land-cover  classification.  Two  ICESat  tracks  are  aligned  to  the  right  of  their  respective 
CORINE  Land  Cover  (CLC)  2000  database  tracks.  All  tracks  are  colored  by  coverage 
classification:  high  vegetation  (green),  urban  (red),  bare/low  vegetation  (white),  or  water 
(blue).  Overall,  Duong’s  classification  results  matched  well  with  the  CLC2000  data,  with 
general  exceptions  in  three  small  areas  (labelled  one  through  three  with  zoomed  in  views 
to  the  left).  Area  one  represents  bare  ground  misclassified  as  urban;  Area  two  is  water 
misclassified  as  bare  ground;  Area  three  is  urban  misclassified  as  forest.  In  each  area  of 
misclassification,  the  waveform  for  one  return  is  illustrated,  it  corresponds  to  the  location 
marked  by  the  white  arrows.  The  black  waveform  is  the  ICESat  data,  and  the  yellow  its 
Gaussian  decomposition. 
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Figure  24.  Classification  Based  on  Waveform  Gaussian  Decomposition  for  Two 
ICESat  Tracks  Superimposed  on  Landsat-7  Imagery. 

Source:  Duong  (2010,  140). 


Research  studies  which  derive  an  understanding  of  Arctic  and  Antarctic  landmass 
elevations,  sea-ice  thickness,  subglacial  features,  and  climate  changes  are  even  more 
impressive.  Kurtz  et  al.  merged  snow  depth  models  with  ICESat  data  for  the  winter  and 
fall  seasons  from  2003  to  2008  (Kurtz  et  al.  2011).  Their  results  indicated  that  the  warmer 
ocean  is  contributing  an  ever-increasing  effect  on  Artie  land  temperatures  due  to  the 
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shrinking  volume  of  sea-ice,  which  acts  as  an  insulator  between  land  and  sea  (Kurtz  et  al. 
2011).  Sorensen  et  al.  used  an  ICESat  measurement-based  procedure  to  predict  the 
changing  mass  distribution  of  the  Greenland  ice  sheet,  concluding  that  annual  losses 
ranged  from  between  191  ±  23  to  240  ±  28  billion  tons  of  ice  (Sorensen  et  al.  2011). 
More  recently,  Howat  et  al.  combined  ICESat  elevation  data,  Worldview  stereopair 
imagery,  Landsat  imagery,  and  ice-penetrating  radar  measurements  to  expose  the  sudden 
and  unexpected  drainage  of  a  subglacial  lake  in  Greenland  (Howat  et  al.  2015). 

Figure  25  provides  an  elevation  map  of  Antarctica  derived  by  Shuman  from 
ICESat/GLAS  collections  over  approximately  one  month:  3  October  to  8  November  2004 
(Shuman  2016).  The  dark  blue  areas  at  0  m  above  sea  level  are  the  ice  sheets.  The 
maximum  elevation,  about  4  km,  defines  the  majority  of  Antarctica’s  eastern  region.  The 
crisscrossing  lines  show  the  satellite’s  ground  tracks  over  time. 

Figure  25.  Antarctica  Elevation  Model  from  GLAS  Collections.  Source:  Shuman 

(2016). 
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Lastly,  scientific  papers  pertaining  to  vegetation  structures,  biomass,  and  cloud/ 
aerosol  distributions  (ICESat’s  secondary  missions)  are  also  detailed  and  extensive.  For 
example,  Pang  et  al.  discovered  that  GLAS  waveform  data  captured  during  the  summer 
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and  early  fall  seasons  contributed  better  intensities  for  estimating  temperate  forest  heights 
in  northeastern  China  (Pang  et  al.  2008).  Goetz  et  al.  analyzed  ICESat  global-altimetry 
and  land- altimetry  datasets,  in  conjunction  with  Moderate  Resolution  Imaging 
Spectroradiometer  (MODIS)  data  and  information  extracted  from  the  Alaska  Burned 
Area  Database,  to  study  disturbed  forests  in  Alaska  (Goetz  et  al.  2010). 

Figure  26  is  a  vertical  profile  of  wildfire  smoke  clouds  in  central  California.  The 
top  left  window  box  exhibits  the  corresponding  ICESat  track  (green  line),  locations  of 
fires  (red  dots),  and  smoke  direction.  In  the  main  image,  the  x-axis  indicates  latitude  and 
longitude,  and  the  y-axis  is  cloud  height  in  km.  Density  is  color-coded  from  dark  blue 
(ambient  atmosphere)  to  white  (dense  clouds/smoke).  Palm  explains  that  the  densest 
sections  of  wildfire  smoke  are  located  primarily  at  ground  level  and  an  altitude  of  about  4 
km  (Palm  2016). 

Figure  26. 

15 

10 

1 

2 
o> 
o 
X 

5 

0 


The  aforementioned  research  publications  represent  just  a  few  examples  of  the 
vast  literature  that  stemmed  from  NASA’s  ICESat  mission.  Undoubtedly,  datasets  from 
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Vertical  Profile  (Height  and  Density)  of  California  Wildfire  Smoke 
Clouds,  28  October  2003.  Source:  Palm  (2016). 


29.119  30.119  31.120  33.120  34.120  36.120 

Latitude,  Longitude  (North.  West) 


^ -  i  I 

low  high 

Density 


its  follow-on,  ICESat-2,  will  motivate  even  more  cryosphere-related,  atmospheric,  and 
forestry  publications,  including  ones  that  make  critical  comparisons  between  current 
ICESat  and  future  ICESat-2  datasets  for  long-term  trend  analyses. 

10.  Mercury  Surface,  Space  Environment,  Geochemistry,  and  Ranging 
Orbiter 

Abshire  pointed  out  that  “in  March  2011,  the  Mercury  Surface,  Space 
Environment,  Geochemistry,  and  Ranging  Orbiter  (MESSENGER)  became  the  first 
spacecraft  in  history  to  orbit  our  innermost  planet”  (Abshire  2011,  41).  He  added  that  the 
instrument  suite  underneath  the  probe  contained  a  magnetometer;  energetic  particle, 
neutron,  gamma  ray,  x-ray,  atmospheric,  surface,  and  plasma  spectrometers;  a  dual  wide- 
and  narrow-angle  imaging  system;  and  the  Mercury  Laser  Altimeter  (MLA)  (Abshire 
2011).  MESSENGER’S  magnetometer  occupies  the  position  at  the  tip  of  the  extended 
pole  in  Figure  27,  and  the  MLA  sits  about  midway  on  the  underside  of  the  spacecraft’s 
main  body. 

Figure  27.  Artist  Impression  of  MESSENGER  Orbiting  Mercury.  Source: 

Hopkins  University  Applied  Physics  Lab/Carnegie  Institution  (2016). 


Enhanced  color  photograph  overlay  taken  near  Hokusai  Crater. 

Cavanaugh  et  al.  processed  MLA’s  range  calculations  and  boresight  pointing 
angles,  in  conjunction  with  spacecraft  orbital  data,  to  precisely  characterize  Mercury’s 
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surface  topography  (Cavanaugh  et  al.  2007).  Additionally,  Neumann’s  analysis  of  surface 
reflectance  measurements  led  to  the  theory  that  anomalous  bright  and  dark  spots  on 
Mercury’s  northern  facing  slopes  indicated  the  possible  presence  of  subsurface  water-ice 
insulated  by  organic  volatiles,  both  deposited  by  foreign  bodies  after  impact  with  the 
planet  (Neumann  et  al.  2012). 

Figure  28  is  a  MLA-derived  topographic  map  of  Mercury  from  4°  South  to  90° 
North  on  a  cylindrical  projection  with  resolution  at  1  km  per  pixel  (NASA/John  Hopkins 
University/Carnegie  Institute  2016).  On  the  elevation  scale  to  the  right,  dark  blue  denotes 
the  lowest  regions  (-5  to  -3  km)  and  brown  denotes  the  highest  points  (2.5+  km),  also: 
light  blue  (-2  km);  green  (-1  km);  yellow  (-0  km);  orange  (1  km);  and  red  (2  km). 


Figure  28.  Topography  Map  of  Mercury’s  Northern  Hemisphere.  Source:  NASA/ 
Johns  Hopkins  University  /Carnegie  Institution  (2016). 


MLA  had  an  8  Hz,  two-stage,  Nd:YAG  laser  transmitter  that  produced  a  1064  nm 
(infrared)  pulse  with  20  mJ  of  energy  and  a  6  ns  pulse-width;  it  weighed  7.4  kg  and 
required  23  watts  of  operating  power  (Abshire  2011).  A  rugged  space-hardened  LiDAR 
indeed,  the  MLA  achieved  more  than  laser  altimetry  and  mapping  of  Mercury. 

During  a  distant- Earth  flyby  in  2005,  MLA  successfully  completed  a  two-way 
laser  linking  experiment  with  the  NASA  Goddard  Geophysical  and  Astronomical 
Observatory  at  a  range  of  over  23.9  million  km  (Smith  et  al.  2006).  Smith  et  al. 
mentioned  that  MLA  transmitted  laser  pulses  to  Earth  during  three  observation  periods 
divided  across  three  nonconsecutive  days  (2006).  They  added  that  Earth  received  brief 
transmissions  of  16  consecutive  signals  on  24  and  27  May  and  a  transmission  of  about  26 
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consecutive  signals  on  the  31st  (the  last  day  of  testing)  (2006).  Meanwhile,  NASA’s 
Earth-based  laser  transmitted  its  own  pulses  toward  MESSENGER,  which  the  MLA 
detector  received  about  0.35  ms  earlier  than  predicted  throughout  the  entire  testing  period 
(2006).  By  comparison,  transmissions  from  MESSENGER  reached  Earth  approximately 
0.34  ms  early  on  the  24th/27th  and  0.14  ms  early  on  the  31st  (2006).  Overall,  laser- 
acquired  range  measurements  between  the  spacecraft  and  Earth  were  accurate  to  within 
±20  cm  at  almost  24  million  km  (2006).  This  ground  breaking  asynchronous  laser 
ranging  and  communications  experiment  “established  a  distance  record  for  laser 
transmission  and  detection  and  accomplished  a  two-way  laser  link  at  interplanetary 
distance”  (Smith  et  al.  2006,  53). 

Once  at  Mercury,  MESSENGER  maintained  a  highly  elliptical  and  nearly  polar 
orbit  with  a  300  km  periapsis  in  the  northern  hemisphere  and  a  15,000  km  apoapsis  in  the 
southern  hemisphere  (Cavanaugh  et  al.  2007).  Cavanaugh  et  al.  noted  that  this  orbit’s 
purpose  was  to  protect  the  valuable  underside  of  the  probe  from  the  intense  solar 
radiation  (2007).  Moreover,  MESSENGER’S  orbital  parameters  ensured  that  the  planet 
itself  shielded  the  scientific  instruments  when  the  probe’s  topside  solar-screen  faced 
away  from  the  sun  (2007).  Of  course,  the  distant  apoapsis  explains  the  sparse  LiDAR 
measurements  south  of  the  equator,  since  MLA’s  maximum  operating  limits  were  1,200 
km  at  nadir  and  800  km  at  slant  range  (Cavanaugh  et  al.  2007).  MESSENGER  exceeded 
its  expected  mission  duration  by  three  years,  ultimately  crashing  into  Mercury’s  surface 
on  30  April  2015  (NASA/Johns  Hopkins  University/Camegie  Institution  2016). 

11.  Cloud- Aerosol  LiDAR  and  Infrared  Pathfinder  Satellite  Observations 

Capturing  first  light  on  7  June  2006,  the  Cloud-Aerosol  LiDAR  and  Infrared 
Pathfinder  Satellite  Observations  (CALIPSO)  spacecraft  is  NASA  Langley  Research 
Center’s  Earth- observing  and  atmospheric-polarization  LiDAR  platform  (Abshire  2011). 
Abshire  states  that  CALIPSO  exploits  three  co-aligned  scientific  instruments:  a  wide- 
field  camera,  an  imaging  infrared  radiometer,  and  a  polarization  LiDAR  (2011).  He  also 
notes  that  a  star  tracker  assembly  ensures  accurate  spacecraft  and  instrument  pointing 
from  CALIPSO’s  sun- synchronous  98.2°  orbit  at  an  altitude  of  700  km  (Abshire  2011). 
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CALIPSO  currently  belongs  to  the  Afternoon  Constellation  (A-train)  that  consists  of  it 
plus  the  Aqua,  Aura,  and  CloudSat  satellites  (Winker  and  Hunt  2007).  Figure  29  is  a 
rendition  of  CALIPSO  and  a  diagram  of  its  sensor  payload. 


Figure  29.  Artist  Rendition  of  CALIPSO  On-Orbit  and  Diagram  of  Payload 

Optics.  Adapted  from  NASA  (2016a). 


The  LiDAR  is  officially  named  the  Cloud-Aerosol  LiDAR  with  Orthogonal 
Polarization  (CALIOP);  CALIOP  shall  not  be  used  hereafter  in  this  review  to  avoid 
confusion  with  the  spacecraft’s  similar  acronym.  In  short,  this  “first  polarization  LiDAR 
in  orbit  is  a  three-channel  elastic-backscatter  LiDAR  optimized  for  aerosol  and  cloud 
profiling”  (Winker  and  Hunt  2007,  1).  Winker  and  Hunt  discuss  how  CALIPSO  still 
welds  one  remaining  of  two  linear-polarizing  Nd:YAG  laser  transmitters  that  produce  a 
combined  pulse  of  1064  nm  (infrared)  and  532  nm  (green)  wavelengths,  with  a  total 
energy  of  220  mJ  (2007).  Furthermore,  the  system  operates  at  20  Hz  with  pulse- widths  of 
20  ns,  and  ground  footprints  are  70  m  in  diameter  and  spaced  333  m  apart  (2007). 

The  avalanche  photodiode  infrared  detector  disregards  the  polarization  of  the 
1064  nm  return  signal  (Winker  and  Hunt  2007).  Winker  and  Hunt  describe  that  the  532 
nm  detector  employs  two  separate  channels  of  photomultiplier  tubes  to  measure  the 
parallel  and  perpendicular  components  of  polarization  within  the  returned  backscatter 
(2007).  They  go  on  to  say  that  the  returning  wavelengths  are  sampled  at  10  Megahertz 
(MHz)  with  an  initial  resolution  of  15  m,  then  receiver  electronics  throttle  resolution  to 

60  m  for  the  infrared  channel  and  30  m  for  the  parallel  and  perpendicular  green  channels 
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(2007).  To  further  ease  the  downlink,  additional  data  reductions  are  performed  by 
onboard  processors  that  average  the  top  sections  of  the  vertical  LiDAR  profiles  to 
compress  the  data  prior  to  transmission  (Winker  and  Hunt  2007). 

Overall,  CALIPSO  continues  to  provide  researchers  with  highly  versatile  LiDAR 
products.  For  example,  Hu  et  al.  used  derived  sea-surface  laser  backscatter  along  with 
wind-speed  measurements  from  Aqua’s  Advanced  Microwave  Scanning  Radiometer — 
EOS  (AMSR-E)  to  examine  global  wave  slope  variances  and  near-surface  wind  speeds 
(Hu  et  al.  2008).  Battaglia  examined  CALIPSO  wavelength  and  polarization  data  for 
three  areas  of  relatively  shallow  water:  Kure  Atoll,  the  Bahamas,  and  Sequoia  National 
Forest  (Battaglia  2010).  He  recognized  a  strong  correlation  between  water  depth  and 
cross-polarization  of  the  532  nm  return,  enabling  him  to  estimate  water  depth  based  on 
the  identifiable  LiDAR  return  signatures  for  shallow  waters  less  than  20  m  deep 
(Battaglia  2010). 

Thanks  to  Hu,  scientists  are  able  to  capitalize  on  a  new  calibration  technique 
called  Elevation  in  Tail  that  substantially  improves  the  accuracy  of  CALIPSO  laser- 
altimetry  data  (Hu  2007).  Elevation  in  Tail  compares  the  transient  responses  of  the 
photomultiplier  detectors  to  the  waveform  shapes  of  pulses  returned  from  the  surface 
(2007).  Further  processing  can  be  applied  to  improve  the  accuracy  of  elevation  data 
products  from  CALIPSO  or  a  comparable  LiDAR  with  polarization  capabilities  (Hu 
2007). 

Lastly,  Figure  30  shows  a  CALIPSO  vertical  LiDAR  profile  of  wildfire  smoke 
clouds  near  Springerville,  Arizona.  The  plumes  tend  to  maintain  an  approximate  altitude 
of  4.5  km,  similar  to  the  earlier  ICESat  profile  of  California  wildfire  smoke.  The  LiDAR 
profile  is  superimposed  on  a  MODIS  Terra  image  captured  nine  hours  later  by  the  Terra 
satellite’s  imaging  spectroradiometer  (NASA  2016a).  This  is  but  one  sample  of  the 
tremendous  amount  of  useful  atmospheric  LiDAR  data  that  CALIPSO  has  collected  over 
the  course  of  its  mission:  10  years  and  counting. 
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Figure  30.  A  CALIPSO  LiDAR  Vertical  Profile  of  Arizona  Wildfire  Smoke. 

Source:  NASA  (2016a). 


12.  Lunar  Reconnaissance  Orbiter/Lunar  Orbiter  Laser  Altimeter 

On  23  June  2009,  the  Lunar  Reconnaissance  Orbiter  (LRO)  achieved  its  elliptical 
and  polar  orbit  around  the  Moon,  where  the  spacecraft  maintains  periselene  and 
aposelene  at  20  km  and  165  km,  respectively  (NASA  2007a).  In  addition  to  mapping 
lunar  topography  and  establishing  a  selenocentric  coordinate  system,  LRO  generates  day 
and  night  temperature  maps,  captures  high-quality  imagery,  and  measures  the  Moon’s 
reflectivity  coefficient  (albedo)  (Keller,  Chin,  and  Morgan  2007).  Keller,  Chin,  and 
Morgan  comment  that  LRO  focuses  specifically  on  the  polar  regions,  searching  for 
perpetually  dark  areas  which  may  contain  traces  of  water-ice  (2007). 

The  satellite  operates  seven  scientific  instruments,  of  particular  interest  is  the 
Lunar  Orbiter  Laser  Altimeter  (LOLA),  a  space  LiDAR  conducting  a  dual  mission  of 
general  exploration  and  preliminary  research  in  support  of  future  robotic  and  manned 
lunar  landings  (Keller,  Chin,  and  Morgan  2007).  Specific  tasks  include  “measuring 
potential  landing  site  slopes  and  surface  roughness”  (NASA  2007a).  Figure  31  is  a 
rendition  of  LRO  on  orbit  (left)  and  a  LOLA-derived  map  of  surface  slopes  near  the  lunar 
south  pole.  The  slope  map,  with  the  south  pole  in  the  center,  shows  poleward  of  75 
degrees  south  latitude  and  is  color-coded  from  no-slope  (purple)  to  30  degrees  of  slope 
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(red).  Generally,  the  areas  of  steeper  slope  (15°-30°)  define  the  outlines  of  impact  craters 
in  this  region. 


Figure  31.  Artist  Rendition  of  LRO  and  Slope  Map  of  the  Lunar  South  Pole. 

Source:  NASA  Goddard  (2016b). 


By  design,  this  LiDAR  is  somewhat  similar  to  the  Mars  Orbiter  and  Mercury 
Laser  Altimeters,  MOLA  and  MLA  respectively  (Keller,  Chin,  and  Morgan  2007). 
However,  LOLA  does  more  than  uphold  the  trend  of  Q-switched,  diode-pumped, 
Nd:YAG  lasers  in  space-based  LiDAR  operations.  With  diffractive  optical  elements 
producing  five  total  beams,  it  is  “the  first  multi-beam  space  LiDAR”  (Abshire  2011,  38). 
Abshire  remarks  that  each  beam  produces  an  individual  5  m  diameter  footprint,  and  the 
five  footprints  are  spaced  25  m  apart  in  the  shape  of  a  cross  (201 1).  This  unique  feature  is 
how  LOLA  can  measure  surface  slope  and  elevation  as  it  orbits  the  Moon  (2011).  More 
specifically,  the  28  Hz  LOLA  functions  at  1064  nm  (infrared)  with  pulse  energies  of  3  mJ 
and  a  5  ns  pulse- width  (201 1).  As  the  five  returning  laser  pulses  enter  the  14-cm  receiver 
telescope,  an  optical  fiber  array  and  timing  electronics  direct  each  to  an  individual  silicon 
avalanche  photodiode  detector  (Abshire  201 1). 
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Figure  32  provides  a  LOLA-derived  elevation  map  of  the  two  sides  of  the  Moon. 
Elevation  in  km  is  color-coded  based  on  a  mean  lunar  radius  of  1,738  km.  This  image 
distinguishes  well  between  the  generally  flatter  “lowlands”  of  the  nearside  and  the  higher, 
yet  heavily  cratered,  farside. 


Figure  32.  Lunar  Reconnaissance  Orbiter:  Detailed  Topography  of  the  Moon. 

Adapted  from  NASA  Goddard  (2016b). 


Now,  one  could  reason  that  LOLA  has  completed  the  earlier  works  of  the  Apollo 

and  Clementine  laser  altimeters.  However,  surface  topography  is  only  one  capability  of 

this  built-for-space  LiDAR:  “Laser  systems  can  be  used  to  track  and  communicate  with 

spacecraft  in  deep  space  to  achieve  better  performance  at  lower  power  and  with  smaller 

size  apertures  than  conventional  microwave  systems”  (Sun  et  al.  2013,  1865) 

Appropriately,  Sun  et  al.  successfully  demonstrated  the  distant-communications  ability  of 

the  LOLA  receiver  by  conducting  one-way  digitally  modulated  (4096-ary  PPM)  laser 

transmissions  employing  Reed  Solomon  coding  during  some  of  the  transmissions  (2013). 

They  mentioned  that  after  encoding,  the  Next  Generation  Satellite  Laser  Ranging  ground 
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station  transmitted  essentially  error-free  data  at  rates  between  200  and  300  bits  per 
second  (2013).  LRO’s  distance  from  Earth  was  about  380,000  km  during  the  experiment, 
and  the  spacecraft  returned  accurate  time  tags  for  every  laser  pulse  that  it  received  (Sun  et 
al.  2013). 

Additionally,  laser  backscatter  from  LOLA  offers  a  means  to  determine  the 
albedo  of  various  lunar  surfaces  (Lucey  et  al.  2014).  Ligure  33  gives  the  whole-surface 
albedo  of  the  Moon  at  1064  nm  on  a  cylindrical  projection.  According  to  Lucey  et  al.,  the 
majority  of  the  surface  has  an  albedo  of  0.29  (2014).  They  continue  by  informing  us  that 
the  high  of  0.33  corresponds  to  the  farside’s  northern  hemisphere  below  the  pole,  and  the 
0.24  low  aligns  almost  completely  with  the  segment  of  the  nearside  that  directly  faces 
Earth  (Lucey  et  al.  2014).  The  LiDAR  systems  of  this  section  each  represent  a  unique  and 
tremendous  contribution  to  the  sciences  of  laser  detection  and  ranging,  space  systems 
operations,  and  remote  sensing.  Still,  this  review  would  not  be  complete  without  a  look  at 
some  foreign  LiDAR  space  platforms  and  a  glimpse  into  the  future. 

Ligure  33.  The  Geometric  Albedo  of  the  Moon  at  1064  nm. 

Source:  Lucey  et  al.  (2014,  119). 
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c. 


OTHER  PLATFORMS 


1.  l’Atmosphere  Par  LiDAR  Sur  Saliout  (“Space  Station  Atmospheric 
LiDAR”) — French/Russian  LiDAR  on  the  Mir  Space  Station 

L’atmosphere  par  LiDAR  sur  Saliout  (ALISSA)  was  a  French  atmospheric 
LiDAR  system  designed  by  the  Service  d’Aeronomie  and  installed  on  the  Priroda  module 
of  the  Mir  Space  Station  in  1996  (Chanin  and  Hauchecorne  1991).  Figure  34  is  a  diagram 
of  the  Priroda  module  and  its  location  on  the  Mir  station.  ALISSA’ s  primary  mission,  as 
described  by  Chanin  and  Hauchecorne,  was  to  capture  detailed  descriptions  of  the 
vertical  structures  of  clouds,  including  the  uppermost  altitudes  of  cloud  tops  at  a 
resolution  of  approximately  150  m  (1991).  They  added  that  the  LiDAR’s  four  Nd:YAG 
lasers  synchronized  at  50  Hz  to  produce  a  combined  40  mJ,  532  nm  (green)  pulse;  the 
lasers  could  also  operate  individually  at  10  mJ  (1991).  The  receiving  telescope  was  40  cm 
in  diameter,  and  reference  prisms  and  a  secondary  mirror  sat  in  front  of  the  telescope  to 
correct  for  mechanical  errors  during  space  flight  (Chanin  and  Hauchecorne  1991). 

Figure  34.  The  Mir  Space  Station  and  the  Priroda  Module  with  the  ALISSA 

LiDAR.  Adapted  from  Dismukes  (2013). 


The  Russian  Institute  of  Applied  Geophysics  maintained  responsibility  for 
ALISSA  (McCormick  2005).  McCormick  emphasized  that  the  system’s  controllers  only 
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employed  test  mode  during  the  first  two  years,  then  ALISSA  endured  a  prolonged 
interruption  in  operations  due  to  issues  internal  to  the  space  station  (2005).  He  went  on  to 
say  that  in  1999,  ALISSA  operated  for  approximately  two  months  and  collected  25 
additional  laser  backscatter  datasets  that  indicated  the  presence  of  multi-layered  clouds, 
the  boundary  layer,  and  orographic  waves  (2005).  In  conjunction  with  the  earlier  success 
of  LITE,  ALISSA  further  demonstrated  “that  even  a  low  power  LiDAR  carried  on  an 
operational  satellite  could  provide  useful  information  for  cloud  description,  including  the 
description  of  the  aerosol  boundary  layer”  (McCormick  2005,  381). 

2.  Hayabusa  Asteroid  Probe — Japan 

In  2003,  JAXA  launched  the  Hayabusa  Asteroid  Explorer,  code  named  MUSES- 
C,  to  study  and  return  samples  from  asteroid  25143  Itokawa  (Japan  Aerospace 
Exploration  Agency  [JAXA]  2016).  The  mission  ended  successfully  in  2010  when 
Hayabusa  returned  to  Earth,  allowing  Japanese  scientists  to  recover  and  authenticate  the 
Itokawa  samples  (JAXA  2016).  See  Figure  35. 


Figure  35.  Hayabusa  (MUSES-C)  and  Sample  from  the  Itokawa  Asteroid. 

Source:  Japan  Aerospace  Exploration  Agency  [JAXA]  (2016). 


Hayabusa  featured  a  Xenon-ion  engine,  a  near-infrared  spectrometer,  a  sample¬ 
collecting  robot,  and  an  autonomous  navigation  system  (JAXA  2016).  According  to 

JAXA,  this  human-free  navigation  system  incorporated  an  optical  camera  and  a  laser 
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altimeter  (JAXA  2016).  In  addition  to  successfully  navigating  the  spacecraft  in  vicinity  of 
the  asteroid,  Hayabusa’s  LiDAR  and  camera  suite  obtained  the  measurements  necessary 
for  scientists  to  calculate  the  bulk  density  of  Itokawa  (Mukai  et  al.  2007).  The 
conclusions  by  Mukai  et  al.  were  a  mass  of  3.54  x  1010  kg  with  6%  error,  a  density  of  1.9 
g/cm-3  with  9%  error,  and  a  bulk  porosity  of  about  40%  (Mukai  et  al.  2007). 

The  Hayabusa  LiDAR  was  a  Q-switched,  laser-diode,  Nd:YAG  transmitter  that 
operated  at  1064  nm  (infrared)  (JAXA  2016).  During  testing,  the  actual  wavelength 
produced  varied  by  about  0.3  nm  per  degree  of  temperature  change,  so  the  developers 
installed  a  thermostat  and  heater  specifically  for  the  LiDAR’ s  transmitter  (Mukai  et  al. 
2007).  Mukai  et  al.  mentioned  that  frequent  temperature  fluctuations  still  affected  the 
laser  during  spaceflight  which  resulted  in  intensity  variations  among  the  outgoing  pulses 
(2007). 

The  Hayabusa  LiDAR’ s  minimum  range  was  50  m  due  to  a  small  angular  offset 
between  the  laser  transmitter  and  receiver  (Abe  et  al.  2006).  Abe  et  al.  remarked  that  the 
predicted  maximum  range  was  50  km,  and  the  first  mission-recorded  laser  return 
occurred  at  a  spacecraft  distance  of  49  km  from  the  asteroid  (2006).  They  described  how 
a  timing  clock  frequency  of  75  MHz  afforded  approximately  1  m  accuracy  at  the 
minimum  operating  range  of  50  m;  laser  accuracy  was  2  m  at  slightly  over  3  km  (2006). 
Abe  et  al.  admitted  that  Hayabusa’s  LiDAR  detected  1,665,548  of  the  4,107,104 
transmitted  pulses,  about  41%,  during  the  asteroid-orbiting  and  pre-landing  phase  from 
10  September  to  25  November  2005  (2006).  Unfortunately,  spacecraft  pointing  accuracy, 
and  with  it  the  number  of  returned  laser  pulses,  decreased  substantially  after  2  October 
due  to  a  major  attitude  control  failure:  two  of  three  reaction  wheels  (Abe  et  al.  2006). 

D.  THE  WAY  AHEAD  FOR  LiDAR  IN  SPACE 

1.  NASA — Next  Generation  Spacecraft  Landing  Integrated  LiDAR 

Space  Shuttle  precision  docking  using  TriDAR  represented  a  major  leap  forward 
in  automated  rendezvous  technology.  Building  on  this  success,  NASA  continues  the 
development  of  a  multi-LiDAR  spacecraft  landing  system  for  future  space  exploration 
missions  (Amzajerdian  et  al.  2015).  Amzajerdian  et  al.  state  that  the  Autonomous 
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Landing  and  Hazard  Avoidance  Technology  Project  (ALHAT)  endeavors  to  successfully 
integrate  its  chosen  trio:  an  advanced  laser  altimeter,  a  Doppler  LiDAR,  and  a  3D  Flash 
Imaging  LiDAR  (2015).  Once  installed  onboard  future  spacecraft,  these  systems  will 
work  together  to  achieve  “autonomous  precision  landing  and  hazard  avoidance 
capabilities”  through  highly  accurate  and  simultaneous  measurements  of  spacecraft 
altitude,  relative  velocity,  orientation,  and  proximity  to  terrain  and  other  navigational 
hazards  (Amzajerdian  et  al.  2015,  1). 

In  2014,  the  ALHAT  integrated  LiDAR  system  demonstrated  itself  on  the 
underside  of  a  Huey  helicopter  and  on  the  rocket-powered  Morpheus  unmanned  test 
vehicle,  as  shown  in  Figure  36.  The  continuous-wave  Navigation  Doppler  LiDAR  uses 
frequency-modulated  laser  pulses  segmented  into  up-ramp  chirp,  constant  frequency,  and 
down-ramp  chirp  to  determine  accurate  ranges  and  relative  velocities  (Amzajerdian  et  al. 
2015).  This  report  comments  that  the  pre-modulated  signal  is  a  single  1.55  micron 
(infrared)  pulse  that  is  split  after  modulation  and  transmitted  simultaneously  in  three 
separate  directions  (2015).  It  also  notes  that  the  Doppler  velocity  during  craft  movement 
is  determined  through  a  comparison  of  the  up-ramp  and  down-ramp  chirp  frequencies; 
the  constant  (center)  frequency  is  used  to  confirm  the  calculated  Doppler  velocities, 
helping  the  system  disregard  erroneous  measurements  (2015).  Next,  a  field- 
programmable  gate  array  applies  fast  Fourier  transforms  to  the  different  frequency 
components  of  all  three  return  signals  in  order  to  calculate  actual  spacecraft  range  and 
line-of-sight  velocity  (2015).  Additional  statistics  for  the  Doppler  LiDAR  are  as  follows: 
velocity  resolution  of  0.2  cm/sec,  line-of-sight  resolution  of  17  cm,  maximum  range  of 
2,200  m,  and  a  data  rate  of  20  Hz  (Amzajerdian  et  al.  2015). 


47 


Figure  36.  Testing  of  Integrated  LiDAR  Landing  Systems  on  a  Morpheus 
Vehicle.  Source:  Amzajerdian  et  al.  (2015,  9). 


The  long-range  laser  altimeter  will  complement  the  Doppler  LiDAR  at  ranges  less 
than  2  km,  but  it  will  function  exclusively  at  greater  distances  (Amzajerdian  et  al.  2015). 
Amzajerdian  et  al.  describe  that  this  multi-return,  automatic  threshold  detection,  laser 
altimeter  operates  at  1.57  microns  (infrared),  and  it  is  capable  of  providing  0.5  m 
precision  at  altitudes  above  20  km  to  correct  for  inertial  measurement  unit  drift  during  the 
approach  phase  of  a  landing  (2015).  Furthermore,  they  report  that  close-up  range 
resolution  is  on  the  order  of  5  cm,  and  the  data  rate  frequency  is  30  Hz  (2015).  The  laser 
altimeter  supplies  accurate  terrain  contours  that  will  enable  successful  terrain-relative 
navigation  with  no  human  inputs  (Amzajerdian  et  al.  2015). 

Lastly,  the  3D  Flash  Imaging  LiDAR  is  designed  for  object  recognition  of  hazards 
greater  than  29  cm  at  a  maximum  slant  range  of  2  km;  each  return  of  the  1.06  micron 
(infrared)  pulse  is  measured  individually  by  a  Readout  Integrate  Circuit  at  the  detector 
(Amzajerdian  et  al.  2015).  Amzajerdian  et  al.  contend  that  this  is  a  tremendous 
improvement  over  traditional  laser-imaging  object-recognition  technologies  which 
employ  a  single  detector  and  scan  using  a  raster  pattern  (2015).  Here,  the  128  x  128  pixel 
detector  array  achieves  an  imaging  rate  of  30  frames  per  second  (2015).  Additional 
statistics  for  the  3D  Flash  Imaging  LiDAR  are  as  follows:  50  mJ  pulse  energy,  8  ns  pulse- 
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width,  100  mm  receiver  aperture,  1,800  m  range  at  nadir,  and  8  cm  range  precision 
(Amzajerdian  et  al.  2015). 

2.  NASA — Advanced  Topographic  Laser  Altimeter  System  and  Swath 

Imaging  Multi-polarization  Photon-counting  LiDAR 

NASA’s  ICESat-2  will  “continue  polar  ice  topography  measurements  with 
improved  precision  laser-ranging  techniques”  as  a  follow-on  to  the  success  of  ICESat 
(Troupaki  et  al.  2015,  1).  Troupaki  et  al.  explain  that  the  Advanced  Topographic  Laser 
Altimeter  System  (ATLAS)  comprises  ICESat-2’s  entire  scientific  instrument  loadout 
(2015).  They  believe  that  the  micro-pulsed,  multi-beam,  ATLAS  will  deliver  a  significant 
improvement  over  ICESat’s  single -beam  GLAS  LiDAR  (2015).  Troupaki  et  al.  mention 
that  this  LiDAR  will  work  at  532  nm  (green)  and  10  Kilohertz  (KHz),  splitting  each  pulse 
into  six  individual  beams,  as  seen  in  Ligure  37  (2015).  Originating  as  a  16-beam  concept, 
the  ATLAS  design  subsequently  dropped  to  nine  and,  then,  to  its  current  six-beam 
configuration  due  to  changing  mission  requirements  (2015).  The  full-width  half¬ 
maximum  pulse  duration  is  less  than  1.5  ns,  and  full  output  power  exceeds  9  watts 
(Troupaki  et  al.  2015). 

Ligure  37.  Drawing  of  ICESat- 2  and  its  Six-beam  ATLAS  LiDAR.  Source: 

NASA  Goddard  (2016a). 


Beam  pair  spacing  is  3.3  km  (NASA  Goddard  2016a). 
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Fibertek’s  photon-counting  ATLAS  overcomes  the  complications  of  orbital 
altitude  versus  laser  power  through  low-photon  detection  ( 6  beams  =  1/6  of  the  total 
power  each,  about  1.5  watts )  (Troupaki  et  al.  2015).  Troupaki  et  al.  state  that  an  Optical 
Filter  Assembly  is  included  to  mitigate  the  increased  detection  of  solar  noise  inherent  to 
photon-sensitive  LiDAR  systems  (2015).  LRO/LOLA  broke  the  decades-long  tradition  of 
spacecraft  welding  single-beam  laser  altimeters.  ATLAS  follows  in  the  footsteps  of 
LOLA  by  ushering  in  a  new  era  of  space-borne  LiDAR:  a  metamorphosis  toward 
advanced  multi -beam  systems  with  photon-counting  capabilities  (Yu,  Stephen  et  al. 
2010). 

ICESat-2  and  ATLAS  are  part  of  NASA’s  organizational  trend  toward  wide- 
swath  photon-counting  LiDARs  for  space-based  topography  and  atmospheric-sensing 
missions  (Yu,  Stephen  et  al.  2010).  SIMPL  is  a  theoretical  satellite  laser-imaging  concept 
predicated  upon  the  actual  Airborne  Instrument  Development  Incubator  program  named 
Slope  Imaging  Multi-polarization  Photon-counting  LiDAR  (same  acronym).  The  airborne 
SIMPL  program  arose  from  NASA’s  desire  to  “advance  key  technologies  and  mitigate 
some  of  the  risk  associated  with  the  LiDAR  Surface  Topography  (LIST)  mission”  (Yu, 
Krainak  et  al.  2010,  1).  Furthermore,  NASA’s  aerial-SIMPL  analyses  of  forest  canopy 
structure  and  available  biomass  were  successful  due  to  the  high  resolution,  dual 
wavelength,  and  polarization  properties  of  the  airborne  SIMPL  LiDAR  (Yu,  Krainak  et 
al.  2010).  Results  from  the  airborne  SIMPL  experiment  have  generated  a  tremendous 
anticipation  for  future  ICESat-2  data  (Harding,  Dabney,  and  Valett  2011). 

Figure  38  shows  point  clouds  for  the  four  polarization  channels  (Green  Parallel, 
Green  Perpendicular,  NIR  Parallel,  NIR  Perpendicular)  of  SIMPL  beam  three  of  four 
during  a  Learjet  overflight  of  Lake  Erie  in  winter  2009  (NASA  2016b).  Per  NASA,  the 
segment  is  an  approximately  260  m  long  strip  across  an  area  with  heavy,  yet  young,  sea- 
ice  coverage  (2016b).  Additionally,  a  video  frame  composite,  taken  during  flight, 
provides  a  visual  representation  of  the  segment  (2016b).  The  four  general  types  of  sea-ice 
present  in  the  data  are  labeled  by  letters:  (A)  dark  nilas  ice,  (B,  E)  fresh  gray- white  ice, 
(C)  skim  ice  over  water,  and  (D)  fractures/leads  of  open  water  within  the  ice  sheet 
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(2016b).  The  perceived  slope  exhibited  by  the  laser  points  stems  from  uncorrected  roll 
bias  during  flight  and  not  from  the  actual  topography  of  the  study  area  (NASA  2016b). 


Figure  38.  Data  Segment  of  Beam  #3  Collected  During  Airborne  SIMPL  Flight 
Over  Lake  Erie,  25  February  2009.  Source:  NASA  (2016b). 
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Although  it  will  lack  the  polarization  aspect  of  SIMPL,  ICESat-2  is  an  attempt  at 
laser  imaging  an  ever-larger  swath  from  space  that  intends  to  satisfy  at  least  some  of 
LIST’S  mission  requirements  and  NASA’s  expectations: 
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(1)  A  medium  cost  mission  to  be  launched  by  NASA  between  2016-2020; 

(2)  a  single-instrument  payload  carrying  an  imaging  lidar  at  low  Earth 
orbit;  (3)  one-time  global  mapping  of  land,  ice  sheet  and  glacier 
topography  and  vegetation  structure  through  the  duration  of  the  mission; 

(4)  observe  topography  and  vegetation  structure  change  through  time  in 
selected  areas;  and  (5)  achieve  five  meter  horizontal  resolution,  0.1  m 
vertical  precision,  and  absolute  vertical  accuracy  for  ground  surface 
topography  including  where  covered  by  vegetation.  (Yu,  Krainak  et  al. 

2010,  1) 

Figure  39  is  a  diagram  of  the  conceptual  principles  behind  a  space-based  Swath 
Imaging  Multi-polarization  Photon-counting  LiDAR  platform.  Numerous  fiber  lasers  and 
single-photon  sensitive  detectors  are  integrated  with  advanced  timing  electronics  to 
generate  enhanced  measurements  of  surface  topography  and  reflectance  (Harding  et  al. 
2007).  The  theoretical  SIMPL  space-LiDAR  will  simultaneously  transmit  laser 
wavelengths  of  532  nm  (green)  and  1064  nm  (infrared),  and  the  received  backscatter  will 
be  depolarized  into  components  parallel  and  perpendicular  to  the  original  pulse  in  order 
to  determine  surface  type  and  roughness  based  on  photon  scattering  due  to  wavelength 
and  polarization  (Harding  et  al.  2007). 


Figure  39.  Satellite  SIMPF  Push-Broom  Faser  Altimeter  Measurement 
Approach.  Source:  Harding  et  al.  (2007,  1-2). 
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3.  NASA — Global  Ecosystem  Dynamics  Investigation  LiDAR 

In  2018,  NASA  intends  to  end  the  almost  17-year  absence  of  a  space  station- 
based  LiDAR  which  began  in  2001  when  Mir,  with  ALISSA  onboard,  reentered  the 
Earth’s  atmosphere.  The  plan  is  to  install  a  laser  system  that  will  measure  the  biomass  of 
the  Earth’s  forests  from  the  ISS,  as  it  orbits  at  approximately  410  km  with  51.65  degrees 
of  inclination  (NASA  Science  Missions  2015).  NASA’s  Global  Ecosystem  Dynamics 
Investigation  LiDAR  (GEDI)  is  specifically  designed  to  measure  the  heights  of  trees  (to 
within  1  m  accuracy)  and  the  amount  of  carbon  stored  within  them  and  other  vegetation 
(2015).  The  data  will  be  analyzed  to  generate  a  better  understanding  of  available  forest 
biomass  and  habitat  quality  (NASA  Science  Missions  2015). 

Probable  reasons  for  the  chosen  host  platform  are  the  operational  and  size 
requirements  of  the  14-pulse  system,  also  the  majority  of  the  world’s  critical  forests  lie 
between  51°  north  and  south  latitudes.  The  ISS  is  the  ideal,  and  possibly  the  only,  space 
platform  that  can  currently  support  the  GEDI.  According  to  NASA,  the  GEDI’s  mass, 
power  requirement,  and  operating  wavelength  are  230  kg,  516  watts,  and  1064  nm, 
respectively  (2015).  NASA  intends  to  park  the  GEDI  at  the  exposed  facility  of  the 
Japanese  Experiment  Module,  red  circle  in  Figure  40  (NASA  Science  Missions  2015). 
Figure  41  is  a  simple  diagram  of  the  system. 


Figure  40.  ISS  and  Future  Location  of  GEDI.  Source:  Dubayah  (2014). 
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Figure  41.  Diagram  of  the  GEDI.  Source:  Dubayah  (2014). 
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Yearly,  the  GEDI  system  will  collect  over  15  billion  laser  waveform  profiles; 
gridded  estimates  of  biomass  have  a  predicted  resolution  of  500  m  (Dubayah  2014). 
GEDI  data  will  be  analyzed  to  identify  and  describe  the  changes  in  forest  biomass  due  to 
prolonged  stressors,  such  as  climate  change,  drought,  forestry,  and  other  human 
interactions  (NASA  Science  Missions  2015).  Three  laser  transmitters  and  splitting  optics 
will  produce  the  GEDI’s  14  parallel  beams,  and  individual  laser  ground  spots  will 
measure  approximately  25  m  in  diameter,  with  60  m  along-track  and  500  m  across-track 
spacing  (creating  a  60  m  X  6.5  km  swath)  (Dubayah  2014).  This  is  illustrated  in  Figure 
42.  Notice  that  the  middle,  nadir-facing,  laser  will  operate  at  full  power  and  has  a 
narrower  across-track  swath  for  stronger  signal  returns.  Conversely,  the  two  outer  lasers 
trade  power  for  wider  coverage,  so  they  will  collect  more,  lower-amplitude,  waveform 
profiles. 
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Figure  42.  GEDI’s  Three-sectioned,  14-Beam,  Swath.  Source:  Dubayah  (2014). 
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4.  Sigma  Space — Single  Photon-counting  3D  LiDAR 

Sigma  Space  Corporation  is  striving  to  push  the  advantages  and  the  altitude  of 
single  photon  detection  ever  higher.  Two  of  its  recently  developed  3D-imaging  LiDARs 
appear  promising  for  future  space-based  laser  altimetry  and  topographic  applications:  the 
High  Resolution  Quantum  LiDAR  System  (HRQLS)  and  the  High  Altitude  LiDAR 
(HAL).  The  medium- altitude  HRQLS  can  detect  single  photons  at  an  altitude  of  over  5 
km,  and  the  HAL  can  currently  operate  at  slightly  above  8.5  km  (Degnan  et  al.  2013). 
The  idea  behind  both  systems  is  to  “provide  contiguous  few  decimeter  resolution 
topographic  coverage  on  a  single  overflight  at  aircraft  speeds  up  to  about  220  knots  (407 
km/hr)”  (Degnan  et  al.  2013,  3). 

Instead  of  single-  or  even  multi-pulse  lasers,  these  systems  emit  hundreds  of 
pulselets  at  frequencies  up  to  32  KHz  with  a  dedicated  timing  mechanism  for  each 
pulselet  (Degnan  et  al.  2013).  Degnan’ s  pulse- widths  are  in  picoseconds  and  the  imaging 
rate  is  millions  of  3D  pixels  a  second  (2013).  Ligure  43  shows  the  large-area  DEM  from  a 
12-hour  airborne  HRQLS  survey  of  a  rugged  and  densely  forested  1,700  km"  section  of 
Garrett  County,  Maryland.  Degnan  et  al.  disclose  the  survey  parameters:  2  km  flight 
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altitude,  0.81  km  single  flight  line  width,  and  an  average  aircraft  velocity  of  278  km/hr 
(Degnan  et  al.  2013). 

Figure  43.  HRQLS  Airborne  LiDAR  Survey  of  Garrett  County,  Maryland 
Imposed  on  a  Google  Earth  Map.  Source:  Degnan  et  al.  (2013,  4). 


The  pulselet  concept  applied  by  Sigma  Space  LiDAR  technologies  stems  from 
Degnan’ s  first-generation  research  at  the  NASA  Goddard  Space  Flight  Center  just  prior 
to  his  retirement  from  NASA  in  2003  (Degnan  2002).  Degnan  concludes  that  for  a 
particular  telescope  aperture  and  a  set  level  of  transmitter  power,  traditional  high-energy 
low-frequency  laser  pulses  perform  significantly  worse  for  laser-altimetry  sampling  than 
a  series  of  low-energy  high-frequency  pulses  used  in  conjunction  with  advanced  timing 
receivers  and  photon-counting  detectors  (2002).  In  the  same  report,  Degnan  affirms  that 
employing  small  receiver  telescopes  improves  the  system’s  ability  to  distinguish  between 
legitimate  ground  returns  and  noise  from  other  photon  sources,  primarily  the  Sun 
(Degnan  2002).  Under  his  direction  as  Chief  Scientist,  Sigma  Space  continues  to  improve 
on  these  3D  photon-counting  LiDAR  systems. 

All  things  considered,  Sigma  Space’s  next  generation  of  single  photon-counting 
LiDAR  is  designed  with  future  space  operations  in  mind.  With  single  photon  sensitivity, 

“contiguous,  high  resolution  topographic  mapping  and  surveying  on  a  single  overflight 
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becomes  possible  with  very  modest  laser  powers  and  telescope  apertures — even  from 
orbital  altitudes”  (Degnan  et  al.  2013,  1).  According  to  Degnan  et  al.,  these  systems  will 
possess  the  higher  spatial  resolutions  and  wider  swath  coverages  required  to  expedite 
large-scale  extraterrestrial  topography  missions  without  overly  sacrificing  precision  and 
accuracy  (2013).  They  also  offer  the  improved  foliage,  water,  and  atmospheric 
penetration  critical  for  Earth-targeted  forestry,  bathymetry,  climate/pollution  monitoring, 
and  meteorology  operations  from  orbiting  satellites  (Degnan  et  al.  2013). 

5.  Multi- wavelength  LiDAR  for  Terrain  Classification 

Although  multi-wavelength  LiDAR  technology  remains  primarily  restricted  to 
low-flying  airborne  platforms,  it  represents  a  growing  field  of  interest  and  capability.  The 
predominance  of  LiDAR-in-space  has  been  monochromatic,  with  exceptions  mainly  for 
atmospheric  purposes:  Shuttle  LITE,  Phoenix  Mars  Lander,  ICESat,  CALIPSO,  and 
ALISSA.  Spaceborne  spectral  LiDAR  systems  may  one  day  collect  new  information  on 
extraterrestrial  landscapes  that  thus  far  have  only  been  visited  by  single-wavelength 
platforms:  Mercury,  Mars,  the  Moon,  and  near-Earth  asteroids.  Understandably,  the 
ground  work  must  occur  first,  which  means  spectral  classification  studies  on  the  unique 
features  of  our  Earth,  including  manmade  objects  and  vegetation. 

Fernandez-Diaz  et  al.  recently  published  a  paper  titled  “Capability  Assessment 
and  Performance  Metrics  for  the  Titan  Multispectral  Mapping  Lidar.”  The  paper 
discusses  two  years  of  Optech  Titan  testing  and  mapping  campaigns  that  “demonstrate 
capabilities  to  classify  land  covers  in  urban  environments  with  an  accuracy  of  90%,  map 
bathymetry  under  more  than  15  m  of  water,  and  map  thick  vegetation  canopies  at  sub¬ 
meter  vertical  resolutions”  (Fernandez-Diaz  et  al.  2016,  8).  Also,  Morsy  et  al.  have 
obtained  classification  accuracies  of  92.51%  using  three  normalized  difference  feature 
indexes  (NDFI)  for  separate  land/water  and  vegetation/built-up  area  analyses  (Morsy  et 
al.  2016).  Their  NDFIoreen-NiR,  NDFIoreen-MiR.  and  NDFInir-mir  indexes  are  calculated 
from  the  spectral  intensity  values  of  a  Titan  dataset  collected  at  two  locations  in  Ontario, 
Canada  (2016).  Morsy  et  al.  convert  the  3D  point  cloud  into  a  raster  with  grids  of  1  m  cell 
size,  assigning  each  grid  a  mean  intensity  value  before  calculating  the  indexes  (2016). 
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They  conclude  that  NDFIcreen-MiR  is  better  for  distinguishing  vegetation,  whereas 
NDFInir  mir  is  better  for  built-up  (urban)  areas  (2016).  All  three  of  then-  NDFIs  perform 
poorly  (32.19%-40.31%)  at  separating  land  from  water,  unless  the  points  are  restricted 
by  “number  of  returns”  (2016).  Inclusion  of  the  “number  of  returns”  attribute 
significantly  reduces  the  misclassification  of  water  points  as  vegetation,  increasing 
success  to  approximately  97%  (Morsy  et  al.  2016). 

In  his  thesis,  Judson  Thomas  applies  a  new  spectral  LiDAR  technique  for  terrain 
classification  using  the  ENvironment  for  Visualizing  Images’  (ENVI)  “n-Dimensional 
Visualizer”  (N-D  VIS)  tool  (Thomas  2015).  Thomas’  data  are  not  rasterized,  and  he  does 
not  average  the  intensities;  the  individual  spectral  returns  for  the  green,  NIR,  and  MWIR 
channels  are  preserved  for  analysis  (2015).  His  primary  focus  centers  on  manmade 
surfaces  in  a  suburb  of  Toronto,  Canada:  various  pavements,  roofing  materials,  road 
paint,  and  railway  tracks  (Thomas  2015).  This  thesis  will  apply  a  similar  process  in  an 
effort  to  spectrally  subclassify  ground  and  roofing  materials  on  the  campus  of  the  Naval 
Postgraduate  School,  where  accurate  ground  truth  is  readily  accessible. 
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III.  DATASET  AND  PREPARATIONS 


A.  INSTRUMENT 

1.  Optech  Titan  Multispectral  LiDAR 

The  Optech  Titan  multispectral  sensor  collected  the  aerial  LiDAR  data  for  this 
thesis.  Titan  employs  three  independent  lasers  with  a  combined  ground- sampling  rate  of 
approximately  1  MHz  (Teledyne  Optech  2015).  It  expands  on  the  traditional  spatial 
applications  of  LiDAR  by  introducing  multispectral  sensing.  According  to  the 
manufacturer,  the  Titan  system  “offers  improved  performance  for  3D  land  cover 
classification,  vegetation  mapping,  shallow  bathymetry,  and  dense  topography” 
(Teledyne  Optech  2015,  3).  This  thesis  strives  to  support  that  claim,  primarily  for  ground 
types  and  roofing  materials.  Table  3  lists  relevant  information  about  the  sensor,  and 
Figure  44  is  a  drawing  of  the  laser  component  from  the  specifications  brochure. 


Table  3.  Some  Applicable  Titan  Attributes. 

Adapted  from  Teledyne  Optech  (2015). 


Wavelength 

Offset 

Beam  Divergence 

Channel  1 

1550  nm  (SWIR) 

3.5°  Forward 

0.35  mrad 

Channel  2 

1064  nm  (NIR) 

0  °  (Nadir) 

0.35  mrad 

Channel  3 

532  nm  (Vis-Green) 

7  °  Forward 

0.7  mrad 

Figure  44.  Titan  Laser  Drawing.  Source:  Teledyne  Optech  (2015). 
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2. 


Dataset  Collection 


In  May  2016,  the  National  Center  for  Airborne  Laser  Mapping  (NCALM),  based 
at  the  University  of  Houston,  flew  the  Optech  Titan  sensor  over  Monterey,  California. 
NCALM  collected  multi-wavelength  and  waveform  LiDAR  data  for  the  Remote  Sensing 
Center  at  the  Naval  Postgraduate  School  (NPS).  The  complete  spectral  dataset  has  23 
flight  lines,  mostly  north-south  tracks,  spanning  from  the  Monterey  Regional  Airport  to 
the  city  of  Pacific  Grove.  There  are  three  laser  format  (LAS)  files  per  flight  line, 
corresponding  to  the  individual  532  nm,  1064  nm,  and  1550  nm  channels.  For  this  effort, 
the  four  north-south  flight  lines  (12  LAS  files  total)  that  overlap  the  NPS  campus  are 
selected:  FL  10,  FL  12,  FL  13,  and  FL  15,  as  shown  in  Figure  45. 

Figure  45.  Titan  Monterey  Dataset — Four-Selected  Flight  Lines  and  Overlap. 


FL’s  10  and  15  correspond  to  the  central  campus,  and  FL’s  12  and  13  correspond 


to  the  western  and  eastern  areas,  respectively.  As  a  reference,  Figure  45  provides  an 
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overhead  view  of  the  four-selected  flight  lines,  colored  by  relative  height,  and  the  flight 
line  overlap.  The  dataset  corresponds  to  terrain  located  in  Universal  Transverse  Mercator 
(UTM)  zone  ION,  and  it  is  referenced  using  the  North  American  Datum  of  1983 
(NAD83)  and  the  North  American  Vertical  Datum  of  1988  (NAVD  88). 

B.  DATA  PREPARATION 

1.  Flight  Line  Channel  Merging,  Boundary,  and  Noise  Removal 

Figure  57  offers  an  overall  roadmap  of  data  preparation  and  processing,  and  all  of 
Section  B  is  summarized  in  the  left  hand  column.  The  first  preparation  step  uses  a  set  of 
locally  developed  ArcGIS  scripts  (Titan  Toolbox)  to  merge  the  three  channels  of  each 
flight  line  into  one  LAS  file  and,  then,  clip  all  of  the  flight  line  boundaries.  The  LAS 
“point  source  id”  field  retains  the  original  channel/flight  line  information  for  each  point. 
The  Titan_boundary_clip  script  is  essential,  because  the  laser  points  become  overly 
spread  out  toward  the  edges  of  a  flight  line — where  the  laser  optic’s  scan  angle  is 
greatest,  up  to  60°  off-nadir  (Teledyne  Optech  2015).  Leaving  the  edge  points  will  cause 
discrepancies  during  the  upcoming  nearest  neighbor  algorithm.  The  data  loss  due  to 
clipping  is  insignificant,  as  the  edges  occur  within  the  areas  of  flight  line  overlap  and  at 
the  absolute  boundaries  of  the  point  cloud  (See  Figure  46). 

Figure  46.  Example  Titan  Monterey  Flight  Line  After  Boundary  Clip. 
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Figure  46  is  an  example  of  a  flight  line  after  running  the  Titan_boundary_clip 
script.  The  points  with  classification  “0”  (red)  will  be  retained  during  the  upcoming 
LASground  script,  however  the  edge  points  (classification  “  1 2' -green)  will  be  discarded. 
Next,  it  is  necessary  to  remove  noise  points  before  proceeding  to  spatially  classify  the 
point  cloud.  Martin  Isenburg’s  LAStools  software  suite  includes  a  script  that 
automatically  removes  point  cloud  noise.  However,  the  raw  Monterey  dataset  is  relatively 
noise  free.  The  manual  removal  of  noise  points  from  all  four  flight  lines  took  less  than  30 
minutes  total  using  Applied  Imagery’s  Quick  Terrain  Modeler  (QTM)  software. 

2.  Standard  Ground,  Building,  and  Vegetation  Classification — Spatial 
Geometry  Only  (LASclassify) 

The  basic  classification  process  begins  with  establishing  the  ground  level  using 
LASground,  another  LAStools’  script.  The  following  input  parameters  produce  the  best 
ground  result  for  the  NPS  campus:  -city  mode,  -offset  12,  -fine.  “City  mode”  applies  a 
step  height  of  35  m  to  help  distinguish  flat  rooftops  from  the  ground.  “Offset  12”  includes 
all  points  12  cm  above  the  initial  estimate  as  ground.  “Fine  mode”  uses  a  greater  number 
of  points  to  create  the  initial  ground  estimate.  After  LASground,  all  of  the  points  are 
classified  as  either  “class  1-unclassified”  ( red  points)  or  “class  2-ground”  (green  points). 
Figure  47  is  a  screenshot  of  the  results  of  LASground  near  the  quad  and  Herrmann  Hall. 

Figure  47.  LASground  Results  for  a  Small  Section  of  the  NPS  Campus. 


62 


Next,  the  above  ground  level  (AGL)  metric  is  calculated  via  the  LASheight  script 
in  order  to  determine  the  height  above  ground  for  each  laser  point.  AGL  is  a  critical 
metric  that  will  be  relied  on  heavily  by  the  upcoming  spectral  analysis  and  the  ENVI 
classifiers.  After  running  LASheight,  the  point  cloud  is  ready  for  the  LASclassify  script, 
which  will  separate  the  unclassified  points  that  correspond  to  buildings  and  vegetation. 
LASclassify  analyzes  neighboring  points  to  identify  planar  structures  (buildings)  and 
rough  surfaces  (vegetation).  The  result  is  a  point  cloud  now  with  four  basic 
classifications:  red-unclassified  (1),  green-ground  (2),  dark  blue-vegetation  (5),  and 
cyan-buildings  (6).  Note  that  this  autonomous,  spatial-only,  classifier  is  limited  and  does 
result  in  discrepancies.  For  example,  in  Figure  48  the  rough  clay-tiled  roof  and  the  walls 
of  Herrmann  Hall  are  misclassified  as  vegetation  in  some  areas  or  left  unclassified. 


Figure  48.  Raw  Results  of  the  FASclassify  Script. 


Thomas’  statement  is  still  true:  “To  provide  increased  fidelity  to  the  follow  on 
spectral  classification,  the  results  of  automatic  classification  tool,  FASclassify,  require 
manual  correction”  (Thomas,  2015,  25).  Here,  the  primary  focus  of  manual  correction  is 
to  fix  the  portions  of  buildings  previously  misclassified  as  vegetation.  Also,  we  classify 
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as  many  unclassified  (red)  points  as  possible  into  either  the  vegetation  or  building 
categories.  This  will  allow  our  spectral  analysis  to  focus  directly  on  the  subclassification 
of  specific  ground  types  and  manmade  materials,  instead  of  trying  to  separate  out 
misclassified  roofing  materials  from  the  vegetation  class.  No  new  classifications  are 
established  manually,  this  is  for  correction  only.  Figure  49  is  the  same  area  as  Figure  48 
after  manual  correction.  Manual  correction  is  accomplished  in  Quick  Terrain  Modeler  by 
highlighting  a  region  with  the  polygon  selection  tool  and  using  the  “set  classification  in 
area”  option. 


Figure  49.  Corrected  Results  of  the  LASclassify  Script. 


Before  proceeding  to  the  nearest  neighbor  script,  the  analysis  area  is  chipped  out 
from  the  four  full-length  flight  lines.  The  analysis  area  contains  the  NPS  campus,  a 
section  of  Monterey  Beach,  and  a  buffer  zone.  The  buffer  zone  includes  CA  Highway  1 
to  the  south,  the  ocean  to  the  north,  and  the  two  side-streets  to  the  east  and  west  of  the 
campus.  Figure  50  is  a  Google  Earth  view  of  the  analysis  area  outlined  in  red.  And, 
Figure  5 1  displays  the  analogous  point  cloud  in  QTM,  colored  by  classification.  Lastly, 
Table  4  break  downs  the  number  of  points  per  classification  within  the  analysis  area. 
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Figure  50.  NPS  Subset  in  Google  Earth.  Adapted  from  Google  Earth  (2016). 


Figure  51.  NPS  Subset  Colored  by  Classification — Post  LASclassify. 
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Table  4.  Point  Classification  Breakdown — Post  LASclassify. 


Classification 

Color 

No.  of  Points 

%  of  Total 

1 — Unclassified 

Red 

3,401,822 

22.64 

2 — Ground 

Green 

6,826,990 

45.45 

5 — Vegetation 

Dark  Blue 

3,970,153 

26.43 

6 — Buildings 

Cyan 

823,837 

5.48 

15,022,802 

100 

Overall  accuracy  is  77.36%  (the  total  percentage  of  points  that  are  not  unclassified). 


3.  Nearest  Neighbor — RGB  and  Conversion  from  LAS  to  ASCII  Files 

A  local  Matrix  Laboratory  (MATLAB)  script  employs  the  OpenTSTOOL’s 
nearest  neighbor  algorithm  for  the  process  of  spectral  incorporation.  This  is  the  principle 
step  of  data  processing,  and  it  makes  future  spectral  analysis  possible.  First,  each  flight 
line  file  is  divided  back  into  three  separate  LAS  files,  one  for  each  channel:  Cl-1550  nm, 
C2-1064  nm,  and  C3-532  nm.  Then,  for  all  of  the  laser  points  in  each  channel,  the 
neighboring  intensities  from  the  other  two  channels  are  located.  The  two  results  are  saved 
in  the  LAS  “Red-Green-Blue”  (RGB)  attributes  of  the  original  point  according  to  the 
following  assignments:  The  nearest  1550  nm  intensity  is  stored  as  the  “Blue”  field.  The 
nearest  1064  nm  intensity  is  stored  as  the  “Red”  field.  And,  the  nearest  532  nm  intensity 
is  stored  as  the  “Green”  field.  Of  course,  each  point  duplicates  its  own  intensity  into  the 
appropriate  RGB  field  before  the  remaining  two  fields  are  populated  by  the  intensity 
values  from  the  other  wavelengths.  Figure  52  illustrates  this  process. 

With  a  spatial  point  density  of  approximately  10-15  points  per  square  meter,  the 
spectral  incorporation  via  nearest  neighbor  is  occurring  with  nearly  coincident  points 
(spatially  and  temporally),  so  discrepancies  from  this  process  are  negligible. 
Additionally,  the  nearest  neighbor  injection  of  spectral  returns  into  the  LAS  “RGB” 
attributes  produces  a  false-color  point  cloud  that  can  be  rasterized,  if  desired,  into  a 
multiband  image.  Figure  53  is  this  point  cloud  in  QTM  colored  by  “RGB.”  Roads, 
buildings,  and  other  manmade  objects  appear  close  to  true  color,  and  vegetation  appears 
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in  varying  shades  of  red  and  pink  (higher  1064  nm  returns).  Areas  of  water,  such  as  the 
ocean,  appear  mostly  dark  due  to  weak  signal  returns  in  all  channels.  The  coastal  surf  and 
the  shallow  NPS  reflecting  pool  appear  green,  indicating  stronger  532  nm  returns. 


Figure  52.  Nearest  Neighbor  Spectral  Incorporation.  Source:  Miller  et  al.  (2016). 


“For  every  point  in  each  channel,  find  the  nearest  neighbor  intensities  from  the  other  2  channels 
and  store  in  the  RGB  fields”  (Miller  et  al.  2016,  14). 


Figure  53.  NPS  Point  Cloud  in  False  Color  RGB — Post  Nearest  Neighbor. 
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4.  Final  ASCII  Preparations  and  Calculation  of  Vegetation  Indices 

In  addition  to  LAS  files,  the  nearest  neighbor  script  outputs  the  point  cloud  as  text 
files,  one  per  flight  line.  These  ASCII  files  provide  the  input  to  the  final  script  of  the  data 
preparation  process.  This  is  necessary  for  two  reasons:  1).  The  ASCII_prep_pro  script 
must  append  several  more  attributes  for  each  data  point,  however  the  LAS  file  structure 
cannot  accommodate  the  additional  fields.  2)  There  is  no  direct  conversion  from  LAS  to 
ENVI  standard  format;  LAS  to  ASCII  to  ENVI  format  is  the  only  way  to  get  the  data  into 
ENVI’s  N-D  VIS  tool. 

ASCII_prep_pro  is  a  locally  developed  Interactive  Data  Language  (IDL)  program 
designed  to  merge  the  separate  ASCII  flight  line  files  into  a  collective  array  of  data 
points,  complete  several  additional  processing  steps,  and  output  the  result  as  a  single 
ENVI  file  set  (.dat,  .hdr,  and  .roi  files).  The  intermediate  ASCII_prep_pro  steps  are  as 
follows:  flag  a  5%  subset  of  points,  calculate  vegetation  indices,  and  append  a  manual 
classification  attribute.  A  reduction  factor  of  20  is  applied  via  random- sampling  code  to 
flag  a  5%  subset  of  the  data  points  evenly  throughout  the  entire  scene.  The  remaining 
95%  of  the  points  are  maintained,  however  the  reduction  flag  allows  them  to  be  separated 
out  during  spectral  analysis  in  N-D  VIS.  This  permits  the  analyst  to  spectrally  classify 
only  the  5%  subset  and,  then,  test  the  abilities  of  various  ENVI  classifiers  on  the 
remaining  95%  (e.g.,  K-means,  Maximum  Likelihood,  or  Spectral  Angle  Mapper).  The 
reduction  flag  is  appended  to  the  data  array  as  attribute  #21  for  each  point:  A  value  of  “0” 
means  that  point  is  part  of  the  95%.  A  value  of  “1”  means  part  of  the  5%  training  subset. 
Figure  54  indicates  how  the  data  points  will  be  separated  by  visualizing  attribute  #21  in 
N-D  VIS;  the  upper  line  is  comprised  of  all  points  with  a  reduction  flag  value  of  “1,”  and 
the  lower  line  contains  all  points  with  a  value  of  “0.” 

Next,  ASCII_prep_pro  calculates  three  vegetation  indices  and  appends  them  to 
the  data  array  as  attributes  #18,  #19,  and  #20.  The  vegetation  indices  will  assist  us  in 
differentiating  points  as  either  vegetation  or  non-vegetation  in  areas  where  it  is  otherwise 
unclear.  Since  Titan  does  not  carry  a  visible-red  laser,  the  three  available  vegetation 
indices  (Table  5)  are  the  Green  Normalized  Difference  Vegetation  Index  (GNDVI), 

Green  Difference  Vegetation  Index  (GDVI),  and  Green  Ratio  Vegetation  Index  (GRVI). 
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Table  5.  Vegetation  Indices.  Adapted  from  Richards  (2013). 


Green  Normalized  Difference 
Vegetation  Index 

GNDVI 

=  (NIR  -  Green)/(NIR  +  Green) 

Green  Difference  Vegetation  Index 

GDVI 

=  NIR  -  Green 

Green  Ratio  Vegetation  Index 

GRVI 

=  NIR/Green 

All  three  are  calculated  for  each  point:  Green  is  the  532nm  intensity  value — stored  in  the  “G” 
field  of  a  point’s  RGB.  NIR  is  the  1064  nrn  intensity  value — stored  in  the  “R”  field. 


Appending  the  manual  classification  attribute  is  an  optional  step.  In  order  to  make 
use  of  this  attribute,  one  must  first  export  the  LAS  files  from  QTM  as  “XYZ”  ASCII 
files.  The  ASCII_prep_pro  script  prompts  for  these  input  “XYZ”  files  after  it  asks  for  the 
post-nearest  neighbor  files.  The  classification  of  each  point  in  the  “XYZ”  files  is  copied 
over  to  the  matching  point  in  the  data  array  and  saved  as  attribute  #17.  This  step 
contributes  a  potential  method  of  ground-truth  comparison,  if  one  is  willing  to  manually 
classify  the  entire  point  cloud.  Again,  we  added  no  new  classifications  to  the  base  four 
during  the  earlier  manual  correction,  making  attribute  #17  not  applicable  to  this  work. 

Lastly,  ASCII_prep_pro  formats  the  ASCII  data  into  an  array  of  the  following 
dimensions:  (1  x  #  of  points  x  #  of  attributes).  The  array  for  the  NPS  analysis  area  is  (1  x 
15,022,802  points  x  21  attributes).  The  21  attributes  for  each  point  are  outlined  in  Table 
6.  The  array  is  automatically  exported  as  an  ENVI  .dat  file  with  accompanying  header 
(.hdr)  and  region  of  interest  (ROI)  files. 


Table  6.  ASCII  Attributes.  Adapted  from  Miller  et  al.  (2016). 


Field  # 

Attribute 

Field  # 

Attribute 

1 

X  (UTM — meters) 

12 

Point  Source  ID 

2 

Y  (UTM — meters) 

13 

Time  (GPS  Tag) 

3 

Z  (meters) 

14 

Red 

4 

Original  Intensity 

15 

Green 

5 

Return  Number 

16 

Blue 

6 

Number  of  Returns 

17 

Manual  Classification 

7 

Scan  Direction 

18 

GNDVI 

8 

Edge  of  Flight  Line 

19 

GDVI 

9 

Classification  (Initial) 

20 

GRVI 

10 

Scan  Angle  Rank 

21 

Reduction  Flag 

11 

Height  (AGL) — 8  bit  relative  scale  from 1 

0  (ground)  to  255 
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5.  Input  to  ENVI’s  N-Dimension  Visualizer  (N-D  VIS) 

The  ENVI  files  are  brought  into  the  N-D  VIS  tool  for  spectral  analysis  by 
selecting  the  “visualize  with  new  data”  option.  N-D  VIS  allows  an  “n”  number  of 
dimensions  to  be  visualized  at  one  time  in  the  display  window.  N-D  VIS  is  designed  for 
traditional  multi-  and  hyperspectral  imagery,  where  an  image’s  wavelength  bands  are  the 
dimensions  for  visualizing  its  pixels.  Here,  we  have  adapted  N-D  VIS  to  our  LiDAR  data 
currently  in  ENVI  file  format,  where  we  will  instead  represent  the  21  attributes  as  the 
dimensions  and  visualize  laser  points  in  this  n-dimensional  space.  Each  laser  point  is 
treated  as  an  “image  pixel”  in  our  data  array  of  15,022,802  lines  (“rows” — one  per  point) 
and  21  samples  (“the  attribute  columns”  of  Table  6): 

Use  the  n-D  Visualizer  to  locate,  identify,  and  cluster  the  purest  pixels  and 
the  most  extreme  spectral  responses  (endmembers)  in  a  dataset  in  n- 
dimensional  space.  The  n-D  Visualizer  was  designed  to  help  you  visualize 
the  shape  of  a  data  cloud  that  results  from  plotting  image  data  in  spectral 
space  (with  image  bands  as  plot  axes). ..When  using  the  n-D  Visualizer, 
you  can  interactively  rotate  data  in  n-D  space,  select  groups  of  pixels  into 
classes,  and  collapse  classes  to  make  additional  class  selections  easier. 

You  can  export  the  selected  classes  to  ROIs  and  use  them  as  input  into 
classification  techniques.  (Harris  Geospatial  Solutions  2017c) 

The  “n-D  Controls”  window  gives  us  the  ability  to  select  which  of  the  21 
attributes  from  Table  6  to  study  at  any  one  time.  Hereafter,  we  generally  use  the  term 
bands  instead  of  attributes  to  maintain  consistency  with  ENVI  and  N-D  VIS  terminology. 
The  goal  of  our  spectral  analysis,  in  Chapter  IV,  is  to  view  the  RGB  spectral  values 
(Bands  14,  15,  and  16)  of  all  points  to  identify  and  reclassify  point  clusters  that 
correspond  to  specific  ground  types,  vegetation,  and  manmade  impervious  surfaces.  This 
represents  a  unique  capability  that  is  not  possible  with  geometric  classifiers,  such  as 
LASclassify.  Additionally,  spatial  attributes  will  be  exploited  as  necessary  to  assist  the 
spectral  classification  process.  For  example,  the  AGL  “height”  (Band  11)  will  be  utilized 
to  separate  the  ocean’s  surface  and  the  seafloor  and  to  distinguish  between  ground  and 
roofing  concretes. 

Two  steps  remain  before  the  spectral  subclassification  can  begin.  First,  data 
reduction  is  accomplished  by  selecting  X  (Band  1)  and  Reduction  Flag  (Band  21)  in  the 
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“n-D  Controls”  window.  In  Figure  54,  the  random  subset  of  points  appears  as  the  upper 
line  (Band  21  value  =  “1”)  and  is  selected  by  circling  it  with  the  cursor.  ENVI  Classifiers 
will  exploit  the  user-classified  5%  subset  as  a  training  region  for  classifying  the 
remaining  95%  of  points  (The  lower  line  in  Figure  54:  Band  21  value  =  “0”). 


Figure  54.  N-D  VIS  Data  Reduction  Via  Reduction  Flag  (Band  21). 


X  (Band  1)  and  Classification  (Band  9)  are  selected  to  assign  colors  in  accordance 
with  the  original  classification  outputs  from  the  LASclassify  script.  The  red,  green,  dark 
blue,  and  cyan  scheme  is  temporarily  preserved  to  maintain  consistency.  Classification 
colors  are  assigned  by  selecting  the  appropriate  color  in  the  “n-D  Class  Controls” 
window  and  circling  the  correct  line  of  points  according  to  Figure  55.  Remember, 
LASclassify  assigns  its  output  classes  a  number  designator  (See  Table  4).  Here,  our 
choice  of  Band  1  and  Band  9  presents  the  points  of  each  numbered  class  as  lines  in 
ascending  order:  unclassified  (0),  ground  (2),  vegetation  (5),  and  buildings  (6).  If  the  user 
has  altered  any  of  the  LASclassify  number  designators  or  manually  created  new  classes 
with  a  different  numbering  scheme,  N-D  VIS  would  still  display  those  classes  in  their 
numerically  ascending  order.  It  is  important  to  note  the  minimum  number  of  dimensions 
that  N-D  VIS  can  display  is  two,  which  is  why  X  (Band  1)  is  also  selected  in  Ligures  54 
and  55.  Again,  the  lines  are  actually  groups  of  points  separated  by  Reduction  Llag  (Band 

21)  in  Ligure  54  and  by  Classification  (Band  9)  in  Ligure  55. 
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Figure  55.  N-D  VIS  Point  Clusters  by  Initial  Spatial  Classification  (Band  9). 
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After  applying  the  classification  colors,  the  elements  of  the  point  cloud  become 
distinguishable  in  XYZ  space,  band  selections:  X  (Band  1),  Y  (Band  2),  and  Z  (Band  3). 
XYZ  space  is  now  the  base  state  for  viewing  the  entire  point  cloud  spatially  in  N-D  VIS. 
Figure  56  shows  the  same  data  points  as  seen  in  Figures  54  and  55,  now  displayed  in  their 
normal  x,  y,  and  z  positions  instead  of  by  classification,  reduction  flag  value,  or  any  of 
the  other  bands  (dimensions).  It  provides  a  recognizable  overhead  view  of  the  campus, 
similar  to  Figure  51.  At  this  point,  the  N-D  VIS  session  is  saved  as  a  .ndv  file  to  establish 
a  quick  restoration  point.  This  is  the  best  way  to  save  and  restore  N-D  VIS  work-in- 
progress,  so  that  it  does  not  have  to  be  exported  back  into  ENVI  or  ASCII  files  until  after 
we  have  run  the  ENVI  classifiers  in  Chapter  V. 
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Figure  56.  NPS  Campus  in  N-D  VIS  XYZ  Space. 
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IV.  PROCESSING  APPROACH 


A.  ANALYSIS  WORKFLOW 

Analytical  software  tools  for  the  terrain  classification  of  spectral  LiDAR  datasets 
remain  virtually  non-existent.  The  Optech  Titan  multispectral  LiDAR  is  still  a  relatively 
new  and  unique  system.  Accordingly,  this  thesis  continues  the  use  of  ENVI’s  N-D  VIS 
tool  to  identify  and  define  spectral  subclasses  from  three-wavelength  Optech  Titan  point 
clouds.  Chapter  III  details  the  conversion  of  the  four-chosen  flight  lines  from  traditional 
LAS  file  format  to  ENVI  file  format.  Figure  57  is  the  analysis  workflow  that  covers  the 
contents  of  Chapters  III,  IV,  and  V:  data  preparations,  spectral  analysis  and  processing 
approach  (i.e.,  the  definition  of  training  regions),  and  classification  results  using  ENVI 
classifiers. 

There  are  four  principle  software  tools  that  make  this  research  effort  possible.  The 
LAStools  suite  contributes  the  critical  preparation  steps  of  ground  determination,  height 
calculation,  and  initial  spatial  classification.  Quick  Terrain  Modeler  and  Cloud  Compare 
(CC)  provide  LiDAR  point  cloud  viewing  and  modeling  in  3D  plus  the  ability  to  view, 
filter,  and  analyze  various  data  attributes,  such  as  intensity,  RGB,  scan  angle,  GPS  time, 
return  number,  and  point  source  id.  The  ENVI  software  suite  contributes  the  N-D  VIS 
tool  for  spatial  and  spectral  analysis,  as  well  as  the  necessary  built-in  classifiers.  ENVI’s 
IDL  interface  allows  for  the  quick  conversion  between  the  ASCII  and  ENVI  file  formats. 
This  gives  us  the  ability  to  export  the  completely  classified  point  cloud  from  N-D  VIS 
back  to  QTM  or  CC  for  viewing. 
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Figure  57.  Analysis  Workflow  for  Spectral  LiDAR  Data. 
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B. 


SPECTRAL  CLUSTERING  IN  N-D  VIS  USING  A  5%  RANDOM  SUBSET 


1.  Establishing  the  Land,  Water,  and  Bathymetry  Interfaces 

First,  we  establish  an  intuitive  color  scheme  in  N-D  VIS:  Ground  points  are  set  to 
sienna,  buildings  to  yellow,  and  vegetation  to  a  medium  green.  Dark  blue  and  cyan  are 
used  to  represent  water  and  bathymetry  (ocean  bottom  measurements),  respectively.  The 
95%  of  points  not  flagged  as  part  of  the  5%  training  subset  are  considered  unclassified 
and  colored  white.  The  unclassified  points  are  hidden  from  view  for  the  duration  of 
analysis  by  unchecking  the  white  class  in  the  “n-D  VIS  Class  Controls”  window.  This 
ensures  that  the  upcoming  spectral-clustering  process  is  not  obstructed  by  the 
disproportionately  large  number  of  unclassified  points  (14+  million). 

Identifying  spectral  clusters  (subclasses)  that  correspond  to  specific  types  of 
ground,  building  materials,  and  vegetation  takes  place  primarily  in  RGB  space:  Red 
(Band  14),  Green  (Band  15),  and  Blue  (Band  16).  Nevertheless,  we  distinguish  the  water 
and  bathymetry  subclasses  using  Z  (Band  3),  original  Intensity  (Band  4),  and  AGL  (Band 
11).  Figure  58  is  a  N-D  VIS  screenshot  showing  the  distinctive  land,  water,  and 
bathymetry  interfaces.  Here,  the  wide  range  of  intensity  values  exhibited  by  the  water 
points  aids  our  subclassification  process.  The  two  upward  spikes  on  the  intensity  axis 
correspond  to  the  ocean  and  the  lake. 


Figure  58.  Water  and  Bathymetry  Distinguished  by  Z,  Intensity,  and  Height. 
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Alternately,  the  water  and  bathymetry  subclasses  can  be  selected  by  visualizing 
two-purely  spatial  attributes:  Y  (Band  2)  and  Z  (Band  3)  creates  the  pseudo-side  view  of 
the  analysis  area  as  seen  in  Figure  59.  The  two  bodies  of  water  occupy  distinct  and 
narrow  ranges  of  Z  values — the  surface  of  calm  water  is  essentially  flat.  This  is  also 
apparent  in  Figure  58.  The  bathymetry  points  appear  in  Figure  59  as  a  slope  below  the 
ocean  level.  Vegetation  and  buildings  are  temporarily  unchecked  in  the  “Class  Controls” 
window,  so  that  they  are  hidden  to  avoid  confusion. 

Figure  59.  Water  and  Bathymetry  Distinguished  Spatially  in  N-D  VIS  in  YZ 

Coordinates. 


We  switch  back  to  XYZ  space  (Bands  1,  2,  3)  to  view  the  land,  water,  and 
bathymetry  interfaces  before  proceeding  to  subclassify  the  ground,  vegetation,  and 
building  classes.  Figure  60  presents  five  of  our  six  classes  in  XYZ  space.  The 
unclassified  95%  (white)  remains  hidden  from  view.  The  ability  to  turn  classes  on  or  off 
as  needed  is  another  extremely  useful  attribute  of  the  N-D  VIS  tool.  Again,  dark  blue  is 
the  water,  cyan  is  the  bathymetry,  yellow  is  the  buildings  class,  green  is  the  vegetation 
class,  and  sienna  is  the  ground  class.  Since  we  are  working  with  only  the  5%  subset  of 

points,  Figure  60  appears  less  dense  than  the  comparable  Figure  56. 
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Figure  60.  NPS  Campus  in  XYZ  Space,  Now  with  Water  and 

Bathymetry  Classes. 


2.  Spectral  Subclassification — Ground  Points 

Next,  ground  points  (sienna)  are  isolated  in  order  to  identify  terrain  categories  in 
RGB  space,  also  referred  to  as  spectral  space.  Figure  61  shows  the  ground  points  only  in 
XYZ  space,  approximately  341,300  points.  Figure  62  is  the  corresponding  spectral  space. 
The  data  points  are  now  represented  in  a  format  that  resembles  a  3D  scatterplot.  The 
visualizer  can  animate  forward  or  backward  through  the  various  plot  orientations,  or  one 
can  use  the  “3D  Drive  Axis”  option  to  manually  maneuver  the  display.  In  Figure  61,  the 
silhouettes  of  the  Del  Monte  Lake  and  several  campus  buildings  are  labelled  to  provide  a 
sense  of  orientation  about  the  scene. 
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Figure  61.  Ground  Points  Only  in  XYZ  Space,  Prior  to  Subclassification. 


Del  Monte 
Lake 


Libra  r\ 


Herrmann  llall 


“King  Hall 


Figure  62.  Ground  Points  Only  in  RGB  Space,  Prior  to  Subclassification. 
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The  spectral  analysis  proceeds  systematically  through  the  data  in  RGB  space  by 
visualizing  the  points  from  various  perspectives.  First,  the  extreme  outliers  are  removed 
to  ease  access  to  the  bulk  of  the  points;  they  are  intentionally  absent  from  Figure  62. 
Next,  points  which  appear  to  form  clusters  are  selected  and  assigned  to  a  new  class 
(color).  The  screen  is  then  switched  back  to  XYZ  space  to  inspect  and  correct  the 
selection  as  necessary.  This  process  is  applied  iteratively  until  all  of  the  ground  points 
have  been  reviewed  and  subclassified.  Figure  63  displays  the  largest  four  of  the  newly 
defined  ground  subclasses.  The  most  appropriate  colors  available  in  ENVI  are  assigned  to 
the  new  subclasses. 


Figure  63.  Ground  Subclasses  in  N-D  VIS:  Dirt,  Grass,  Concrete,  and  Asphalt. 
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The  grass/ivy  (bright  green)  points  have  strong  spectral  returns  in  the  red  band 
(1064  nm).  Thus,  they  occupy  the  highest  values  on  axis  14.  The  sidewalk  concrete 
(thistle — gray)  is  relatively  bright  and  nearly  even  in  the  red,  green  (532  nm),  and  blue 
(1550  nm)  bands.  Comparatively,  road  asphalt  (thistle  1)  appears  dark  in  all  three  bands. 
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The  remaining  sienna-colored  points  correspond  to  areas  of  dirt,  sand,  and  mulch,  which 
seem  to  be  spectrally  indistinguishable  in  the  three  available  wavelengths.  Additionally, 
two  smaller  ground  subclasses  are  identified  in  Figure  64.  The  turf  (fake  grass)  near  the 
NPS  flagpole  (colored  sea  green)  and  the  encompassing  red-brick  (red  3)  traffic  circle  are 
the  two  smallest  subclasses  distinguishable  in  RGB  space.  In  Figure  64,  the  grass/ivy 
subclass  remains  displayed  as  a  reference,  but  it  is  necessary  to  hide  the  dirt,  concrete, 
and  asphalt  subclasses  to  view  the  turf  and  brick  spectral  clusters. 

Figure  64.  Ground  Subclasses  in  N-D  VIS:  Grass,  Red  Brick,  and  Turf. 
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The  ground  class  is  now  divided  into  six  spectrally  defined  subclasses:  dirt/sand/ 
mulch,  grass/ivy,  sidewalk  concrete,  road  asphalt,  red  brick,  and  turf.  Figure  65  is  a 
spectral  plot  of  these  ground  subclasses  created  using  the  average  (mean)  RGB  values  of 
each  cluster.  The  x-axis  is  comprised  of  the  three  wavelengths:  532  nm,  1064  nm,  and 
1550  nm.  The  y-axis  is  return  intensity  in  arbitrary  units;  intensities  are  represented  as  a 
12-bit  relative  scale  from  0  to  4096  (2  raised  to  the  12th  power).  None  of  the  mean 
intensities  for  the  ground  subclasses  exceed  2,500.  From  the  3-band  spectral  plot,  the  red 
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brick  and  the  dark  road  asphalt  are  the  two  most  similar  spectra.  Also,  note  the  vast 
differences  between  grass  and  green  turf,  especially  in  the  red  (1064  nm)  and  blue  (1550 
nm)  bands.  Concrete  and  road  asphalt  have  similarly  shaped  profiles;  however,  the 
concrete  is  nearly  twice  as  bright  in  all  bands.  This  corresponds  with  the  locations  of  the 
concrete  and  asphalt  clusters  in  RGB  space  (Figure  63). 


Figure  65.  Spectral  Plot  of  Ground  Subclasses  from  Mean  RGB  Intensities. 
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The  AGL  for  all  ground  points  is  0. 


Before  proceeding  to  subclassify  the  vegetation  and  buildings,  the  newly  defined 
ground  types  are  reviewed  altogether  in  XYZ  space.  The  analysis  is  temporarily  paused 
here  to  permit  a  ground  truth  walk- around  of  the  NPS  campus.  This  occurred  in  late 
October  2016,  approximately  five  months  after  data  collection.  There  was  no  noticeable 
change  in  any  of  the  areas  of  concrete,  asphalt,  turf,  or  the  red  brick  traffic  circle.  Due  to 
a  couple  of  recent  rains  and  gardening,  some  of  the  real-world  areas  of  grass  appear 
larger  and  lusher  than  in  the  dataset.  Thus,  Google  Earth  imagery  from  April  2016  (about 
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one  month  before  the  LiDAR  collect)  is  used  as  a  backup  means  of  ground-truth 
comparison  for  grass  versus  dirt  areas.  Refer  to  Figure  50. 

Figure  66  displays  the  six  ground  subclasses  in  XYZ  space.  Again,  several 
landmark  features  and  building  silhouettes  are  labelled  as  a  guide.  Since  N-D  VIS  only 
offers  three  shades  of  thistle  (gray),  the  boundaries  between  sidewalk  concrete  and 
asphalt  are  difficult  to  discern  in  the  image.  During  analysis,  this  is  overcome  by 
selecting  only  one  class  to  view  at  a  time  by  means  of  the  control  panel.  Lastly,  it  is 
readily  apparent  that  the  grass  on  the  softball  fields  and  around  the  front  of  Herrmann 
Hall  is  the  healthiest  and  densest  on  campus. 


Figure  66.  Ground  Points  Only  in  XYZ  Space,  Post-Spectral  Subclassification. 
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3.  Spectral  Subclassification — Vegetation  Points 

The  primary  intent  of  this  thesis  is  to  spectrally  distinguish  specific  ground  and 
building  materials.  Therefore,  the  vegetation  class  received  the  least  amount  of  analysis. 
A  brief  look  at  the  vegetation  points  in  RGB  space  reveals  two  large  clusters.  The 
spectral  cluster  with  the  greater  intensities  in  all  bands  roughly  corresponds  to  shrubs, 
shorter  cacti,  and  succulents.  The  other  cluster,  with  lower  RGB  intensities,  represents 
taller  bushes  and  trees. 

Figure  67  displays  the  two  newly  defined  vegetation  subclasses  in  RGB  space. 
Shrubs  are  colored  using  “green  2”  and  trees  are  the  slightly  darker  “green  3.”  Although 
grass  technically  belongs  to  the  ground  subclasses,  it  is  included  in  Figure  67  to  show  that 
it  has  stronger  returns  in  the  red  band  (1064  nm)  than  the  taller  vegetation.  The  red  lines 
are  dividers  among  the  three  classes,  not  axes. 

Figure  67.  Vegetation  Subclasses  Compared  to  Grass  in  N-D  VIS. 


We  failed  to  discover  any  distinctive  clusters  corresponding  to  specific  species  of 


trees  or  shrubs.  There  does  not  appear  to  be  enough  information  with  only  three  spectral 

bands  to  distinguish  between  deciduous  or  evergreen  plants,  so  separating  vegetation  at 
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the  species  level  seems  even  less  likely.  For  this  reason,  a  spectra  plot  is  not  included  for 
the  vegetation  subclasses,  and  the  analysis  moves  forward  to  the  last  category:  building 
materials  and  powerlines. 

Figure  68  highlights  the  two  new  vegetation  subclasses  alone  in  XYZ  space. 
Trees  are  the  darker  green  and  generally  appear  denser  than  the  shrubbery.  The  grass 
subclass  is  now  hidden  to  prevent  confusion,  because  the  shades  of  green  in  N-D  VIS  are 
quite  similar.  Several  noteworthy  features  are  labelled,  such  as  the  NPS  Arizona  Cactus 
Garden  and  the  hedges  surrounding  the  flagpole  in  front  of  Herrmann  Hall. 

Figure  68.  Vegetation  Points  Only  in  XYZ  Space, 

Post-Spectral  Subclassification. 
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Mean  AGL  for  shrubs  is  16.5.  Mean  AGL  for  trees  is  127.6. 
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4.  Spectral  Subclassification — Building  Points 

Lastly,  this  analysis  will  apply  the  N-D  VIS  spectral-clustering  technique  to  the 
campus  buildings  and  powerlines.  Figure  69  exhibits  the  starting  point  for  buildings  in 
XYZ  space  (left)  and  the  corresponding  RGB  space  (right).  Several  features  are  labelled 
in  the  XYZ  image,  including  the  powerlines  that  run  along  Del  Monte  Avenue  just 
outside  of  the  campus  fence.  Clusters  of  points  are  more  apparent  in  the  starting  RGB 
image  for  buildings  than  they  were  for  the  ground  and  vegetation  classes.  This  proved 
beneficial,  because  the  subclassification  process  for  buildings  took  approximately  half  the 
time  as  compared  to  the  work  on  the  ground  points.  Additionally,  buildings  are 
successfully  subclassified  into  eight  categories,  as  compared  to  the  six  for  ground  and 
two  for  vegetation.  Building  points  that  cannot  be  spectrally  identified  as  a  specific  type 
of  material  are  retained  in  the  general  building  (yellow)  subclass. 


Figure  69.  NPS  Buildings  in  XYZ  and  RGB  Space, 
Prior  to  Spectral  Subclassification. 


The  first  building  materials  that  we  identify  from  Figure  69  are  two  types  of 
concrete:  light  and  dark.  Light  concrete  (thistle  2)  is  the  brightest  material  in  all  three 
spectral  bands,  and  it  appears  to  also  include  white-plaster  surfaces.  Dark  concrete 
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(yellow  3 — “gold”)  has  a  similar  RGB  proportionality  to  light  concrete,  but  it  maintains 
lower  intensities  across  the  board.  See  Figure  71. 

Clay-tile  roofing  shingles  (coral — color,  not  the  organism)  are  approximately 
twice  as  bright  in  the  red  (1064  nm — NIR)  and  blue  (1550  nm — SWIR)  bands  as  they  are 
in  the  green  (532  nm — visible  green).  Accordingly,  they  are  almost  as  readily  identifiable 
as  the  light  concrete.  Once  the  clay  shingles  subclass  is  separated  out,  it  is  possible  to 
distinguish  another  type  of  roofing  material  that  is  spectrally  similar:  red  lava  rocks  and 
red  basalt-felt  shingles  (orange  3).  This  cluster  sits  just  below  the  clay  shingles  in  RGB 
space,  because  it  is  slightly  dimmer  in  all  three  bands.  Red  lava  rocks  and  red  basalt-felt 
composite  shingles  are  incorporated  into  one  class,  since  they  are  similar  in  both  color 
and  volcanic  origin. 

Three  final  material  clusters  are  identified  with  relatively  low  RGB  intensities. 
These  are  “darker”  than  all  of  the  aforementioned  building  subclasses.  Tan/gray  roofing 
shingles  (maroon  2 — ’’magenta”),  dark  asphalt  shingles  (purple  3),  and  powerlines 
(orchid — ’’light  pink”)  represent  the  least-sampled  building  subclasses  in  the  dataset. 

Figure  70  provides  two  N-D  VIS  screenshots  of  the  new  building  spectral 
subclasses  in  RGB  space.  In  the  right  image,  the  dark  concrete,  clay  shingles,  and  tan / 
gray  shingles  subclasses  have  been  disabled  in  order  to  expose  the  smaller  clusters  for 
dark  asphalt  shingles  and  red  basalt/lava  rock.  Figure  71  is  the  corresponding  spectral 
plot  with  the  three  wavelengths  on  the  x-axis  and  return  intensities  on  the  y-axis.  Y-axis 
values  are  the  mean  RGB  intensities,  in  arbitrary  units,  for  each  subclass.  The  general 
building  (yellow)  subclass  is  included  in  the  plot,  but  would  be  distracting  in  Figure  70. 
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Figure  70.  New  Building  Clusters  in  N-D  VIS,  Less  the 
General  Building  Subclass. 
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Figure  71.  Spectral  Plot  of  Building  Subclasses  from  Mean  RGB 

Intensities  and  Average  AGLs. 
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Figure  72  presents  the  new  building  subclasses  in  XYZ  space.  The  primary 
members  of  the  new  light  concrete  (thistle  2)  subclass  are  King  Hall,  Glasgow  East, 
Halligan  Hall,  most  of  Watkins  Hall,  some  of  the  library  roof,  and  the  plastered  walls  of 
Herrmann  Hall.  Darker  concrete  (yellow  3)  occupies  the  rooftops  of  Bullard  Hall, 
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Spanagel  Hall,  main  Glasgow  Hall,  and  the  flat  section  of  the  Naval  Exchange  (NEX). 
Clay-tile  shingles  (coral)  top  Herrmann  Hall,  the  adjacent  small  houses,  the  peaked 
section  of  the  gym  roof,  and  the  guard  shacks  at  the  two  campus  entrances. 

Red  basalt/lava  rock  (orange  3)  roof  coverings  are  present  on  Ingersoll  Hall,  part 
of  the  library,  the  flat  section  of  the  gym,  and  the  U.S.  Post  Office.  Tan/gray  shingles 
(maroon  2)  can  be  found  on  a  structure  attached  to  the  backside  of  Herrmann  Hall  near 
the  loading  dock  and  in  a  neighborhood  to  the  east  of  campus.  Dark  asphalt  shingles 
(purple  3)  are  present  on  one  house  (#283)  near  the  library  and  Root  Hall  and  in  the  same 
neighborhood  to  the  east. 

Figure  72.  Building  Points  Only  in  XYZ  Space,  Post-Spectral  Subclassification. 
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Ground  truth  involved  accessing  the  upper  balconies  of  Herrmann  Hall  and  the 
rooftops  of  Spanagel  and  Ingersoll  Halls  with  binoculars  to  determine  the  composition  of 
the  other  rooftops  not  visible  from  the  ground.  Since  the  neighborhood  to  the  east  of 
campus  is  not  within  the  primary  area  of  concern,  Google  Earth  is  used  to  approximate 
the  colors  and  types  of  roofing  shingles.  The  points  remaining  in  the  general  building 
(yellow)  subclass  appear  to  be  mostly  metal,  wood,  and  concrete  outliers  that  do  not  fall 
into  the  dark  or  light  concrete  subclasses.  Metal  surfaces,  such  as  the  rooftop  over  the 
NEX  gas  pumps,  have  a  wide  variance  of  intensity  values  and  do  not  form  distinct 
clusters  in  RGB  space.  This  phenomenon  is  likewise  exhibited  on  a  grander  scale  by  the 
water  points,  and  it  is  why  water  requires  subclassification  using  spatial  attributes  in 
addition  to  intensity  (Figure  58).  In  Figure  73,  the  general  building  subclass  is  hidden 
from  display.  This  allows  for  a  better  view  of  the  primary  material(s)  on  each  roof.  Also, 
Figure  73  exclusively  shows  the  NPS  campus. 
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Figure  73.  Roofing  Materials  of  the  NPS  Campus. 
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5.  Summary  of  New  Spectral  Classifications 

By  combining  spatial  and  spectral  attributes  in  N-D  VIS,  this  analysis  session 
produces  19  pertinent  classes  and  one  discarded  class  of  extreme  spectral  outliers.  The  19 
classes  include  the  bulk  95%  (unclassified  points)  and  the  new  ground,  vegetation,  and 
building  subclasses.  These  classes  are  exported  from  N-D  VIS  as  regions  of  interest 
(ROIs)  and  are  ready  to  be  utilized  as  inputs  for  the  ENVI  classifiers  described  in  Chapter 
V.  Figure  74  is  a  screenshot  of  the  newly  classified  5%  training  subset  in  XYZ  space.  The 
shrubs  subclass  is  included,  but  not  labeled  in  the  key.  This  is  because  the  large  extent  of 
Figure  74  makes  them  difficult  to  discern  from  the  similarly  colored  and  more-important 
grass.  Also,  Figure  74  does  not  display  the  unclassified  95%  or  the  extreme  spectral 
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outliers.  Table  7  provides  summary  information  to  accompany  Figure  74;  it  includes  the 
unclassified  points  and  the  outlier  class  for  a  general  comparison. 


Figure  74.  End  of  Spectral  Analysis  in  N-D  VIS — New  Subclasses 

in  XYZ  Space. 
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Table  7.  List  of  Spectral  Subclasses  Generated 
Using  the  5%  Random  Subset. 


Subclass 

Assigned  Color 

Previous  Class — Pre  N-D  VIS 

Number  of  Points 

Unclassified 

White 

All — Randomly  Sampled 
as  described  in  CH  III 

14,338,344 

Extreme  Outliers 

Orange  1 

All 

4,603 

Water 

Blue 

Unclassified  (1) 

107,059 

Bathymetry 

Cyan 

Unclassified  (1) 

11,750 

Dirt,  Mulch,  &  Sand* 

Sienna 

Ground  (2) 

196,044 

Grass  &  Ivy 

Green 

Ground  (2) 

26,204 

Sidewalk  Concrete 

Thistle 

Ground  (2) 

31,872 

Road  Asphalt 

Thistle  1 

Ground  (2) 

72,356 

Red  Brick 

Red  3 

Ground  (2) 

298 

Turf 

Sea  Green 

Ground  (2) 

122 

Shrubs 

Green  2 

Vegetation  (5) 

5,995 

Trees 

Green  3 

Vegetation  (5) 

190,083 

General  Building 

Yellow 

Building  (6) 

7,025 

Light  Concrete 

Thistle  2 

Building  (6) 

6,835 

Dark  Concrete 

Yellow  3 

Building  (6) 

12,078 

Powerlines 

Orchid 

Unclassified  ( 1)  &  Building  (6) 

615 

Clay  Shingles 

Coral 

Building  (6) 

6,799 

Tan  &  Gray  Shingles 

Maroon  2 

Building  (6) 

1,313 

Dark  Asphalt  Shingles 

Purple  3 

Building  (6) 

819 

Red  Basalt  &  Lava  Rock 

Orange  3 

Building  (6) 

2,588 

*Note:  After  the  completion  of  this  chapter,  we  realized  that  sand  could  be  separated  in  RGB  space  from 
the  larger  dirt,  mulch,  and  sand  subclass.  The  coastal  sand  of  Monterey  Beach  forms  a  spectral  cluster  of 
14,179  points;  this  leaves  the  dirt  and  mulch  subclass  at  181,865  points.  The  sand  subclass  is  not  shown  in 
Figure  74,  however  it  will  be  incorporated  into  the  classifiers  and  displayed  in  Chapter  V. 


94 


V.  CLASSIFICATION  ANALYSIS  AND  RESULTS 


A.  GENERAL  PARAMETERS  AND  UNSUCCESSFUL  CLASSIFIERS 

Before  displaying  the  classification  results,  it  is  necessary  to  explain  the 
methodologies  and  input  parameters  applicable  to  all  of  the  ENVI  classifiers 
demonstrated  in  this  chapter.  Additionally,  we  will  briefly  discuss  the  classifiers  that  fail 
to  produce  acceptable  results.  Section  B  continues  with  the  somewhat- successful  K- 
means  unsupervised  classifier,  and  Section  C  describes  the  more-successful  Maximum 
Likelihood  (ML)  supervised  classifier.  Lastly,  Section  D  illustrates  limited  post¬ 
processing  (clean-up)  applied  to  the  ML  results. 

Lirst,  we  consider  an  important  question:  How  do  the  results  differ  if  a  classifier  is 
run  on  just  a  single  flight  line  versus  the  four-combined  flight  lines  of  the  entire  NPS 
campus?  To  answer,  we  always  run  two  concurrent  instances  of  each  classifier  with 
identical  input  parameters — one  instance  classifies  only  a  single  flight  line,  and  the  other 
instance  classifies  the  entire  NPS  campus.  We  choose  flight  line  10  for  all  single  flight 
line  instances,  because  it  offers  the  longest  north-south  run  down  the  center  of  the 
campus.  Plight  lines  12  and  13  are  heavily  cropped  and  exclude  Herrmann  Hall,  and 
flight  line  15  is  essentially  identical  to  10. 

Second,  we  discover  early  in  the  classifier  experiments  that  the  wide  range  of 
intensity  values  exhibited  by  water  points  feeds  excessive  noise  into  the  process.  Water 
points  are  heavily  misclassified  as  vegetation  and,  to  a  lesser  extent,  light  building 
concrete.  To  circumvent  this  problem,  we  need  an  ENVI  mask  file  that  forces  the 
classifiers  to  ignore  the  water  and  bathymetry  points,  both  as  ROIs  (input  classes)  and  as 
points  to  classify.  So,  we  turn  on  the  unclassified  (white)  class  in  N-D  VIS  and  manually 
select  all  water  and  bathymetry  points  via  the  manner  displayed  in  Pigure  58.  Then,  we 
build  a  mask  that  will  exclude  the  water  (blue)  and  bathymetry  (cyan)  classes  by 
choosing  the  “select  areas  off’  option  in  the  ENVI  “Mask  Builder”  tool.  These  masked 
points  will  be  exported  as  a  combined,  yet  otherwise  unaltered,  class  ignored  by  all  of  the 
upcoming  classifiers. 
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Next,  we  experiment  with  various  K-means  and  ML  spectral  input  bands  in 
addition  to  the  height  above  ground  (AGL)  and  Red-Green-Blue  (RGB).  GNDVI  and  the 
other  vegetation  indices  do  not  contribute  a  noticeable  improvement  to  the  results  of  any 
classifier.  Including  the  “number  of  returns”  metric  has  no  impact  on  K-means,  but  it 
does  alter  the  ML  output.  For  ML,  the  “number  of  returns”  metric  improves  the  results 
for  some  classes,  primarily  single-return-only  classes  (e.g.,  rooftop  materials,  turf,  red 
brick,  road  asphalt,  and  sidewalk  concrete).  However,  it  significantly  degrades  the  results 
for  multi-return  classes,  such  as  shrubs  and  trees.  ENVI  classifiers  are  programmed  to 
calculate  statistics  on  continuous  variables  with  expansive  ranges  (e.g.,  AGL  [8-bit  range] 
and  RGB  [12-bit  range]),  but  number  of  returns  is  technically  a  discrete  variable  with  a 
finite  set  of  possible  values:  1,  2,  3,  or  4.  Later,  we  will  make  use  of  the  “number  of 
returns”  metric  with  Quick  Terrain  Modeler’s  “Multivariate  Filtering”  tool,  where  we  can 
selectively  apply  filtering  by  number  of  returns  to  only  the  classes  that  will  benefit — as 
opposed  to  the  classifiers  which  apply  it  indiscriminately.  This  is  the  post-processing 
technique  discussed  in  Section  D,  the  only  post-processing  that  we  perform  on  the  data. 

Lastly,  the  Spectral  Angle  Mapper  (SAM)  and  Spectral  Information  Divergence 
(SID)  supervised  classifiers  produce  extremely  poor  results.  Like  ML,  SAM  and  SID 
generate  an  output  class  for  all  input  classes  from  the  5%  training  subset,  however  there 
is  no  correlation:  All  of  the  SAM/SID  output  classes  contain  points  randomly  dispersed 
throughout  the  scene  that  in  no  way  resemble  the  input  classes.  The  SAM  classifier 
defines  classes  according  to  n-dimensional  spectral  angles,  and  three  wavelengths  does 
not  appear  to  be  enough  for  SAM  to  function  properly.  The  SID  classifier  works  based  on 
n-dimensional  spectral  divergence,  similar  to  SAM’s  spectral  angles.  SID  is  also  hindered 
by  only  three  input  bands.  Moreover,  SAM  and  SID  have  no  capacity  to  consider  the 
AGL  or  any  spatial  attribute;  they  are  strictly  spectral  (RGB  only  for  this  data)  classifiers. 
The  point  clouds  classified  by  SAM  and  SID  are  unsalvageable  and  receive  no  further 
attention. 
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B.  K-MEANS  UNSUPERVISED  CLASSIFICATION 

Since  K- means  is  an  unsupervised  classification  technique,  it  disregards  the  5% 
training  subset.  Instead,  K-means  classifies  the  data  into  its  own  algorithm-defined 
classes  up  to  a  user-specified  maximum.  This  classifier  begins  by  assuming  that  all  class 
means  within  the  data  space  are  distributed  evenly,  then  pixels  are  grouped  into  the  class 
with  the  nearest  mean  via  a  minimum  distance  approach  (Harris  Geospatial  Solutions 
2017a;  Tou  and  Gonzalez  1974).  This  is  an  iterative  process  that  recalculates  class  means 
and  regroups  pixels  until  the  algorithm  reaches  the  user-set  maximum  number  of 
iterations  or  the  overall  change  in  pixel  reclassification  reaches  a  minimum  limit  (default 
setting  is  <  5%)  (Harris  Geospatial  Solutions  2017a). 

We  include  K-means  in  this  work  primarily  to  show  the  limits  of  autonomous 
spectral  classification — similar  to  how  we  described  the  limits  of  autonomous  spatial- 
only  classification  in  Chapter  III.  Although  there  exists  no  quantitative  method  to 
compare  the  unsupervised  results  of  K-means  to  the  supervised  results  of  ML,  it  is  clear 
that  ML  is  the  better  classifier  for  our  spectral  LiDAR  dataset.  This  justifies  the  time  and 
effort  spent  on  the  spectral  training  work  in  Chapter  IV.  Still,  the  K-means  classifier 
performs  better  than  expected  for  both  a  single  flight  line  and  the  combined  flight  lines  of 
the  entire  campus.  The  following  represent  our  best  K-means  results,  with  these 
parameters:  maximum  class  limit  of  nine  (arbitrarily  chosen  as  a  balance  between 
LASclassify’s  four  spatial  and  ML’s  15+  spectral  output  classes);  input  bands  11  (AGL), 
14  (Red),  15  (Green),  and  16  (Blue);  maximum  number  of  iterations  20  (no  improvement 
beyond  this  value);  and  the  default  change  threshold  percentage  of  <  5%. 

Table  8  is  the  point  classification  breakdown  of  K-means  results  for  the  single 
flight  line  and  entire  campus  instances.  It  provides  the  class  names  and  color-codes 
applicable  to  all  upcoming  K-means  figures.  The  two  instances  of  K-means  produce  the 
same  10  classes:  the  nine  defined  by  K-means  and  the  masked  water/bathymetry  class. 
However,  the  individual  point  classifications  differ  between  the  two  K-means  point 
clouds,  especially  at  the  flight  line  edges. 
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Table  8.  K-means  Point  Classification  Breakdown  and  Approximate  Percentages. 


K-means  Class 

(Number  and  Color) 

Name*  and  Description 

Single  Flight  Line 

Entire  Campus 

1 — Dark  Green 

Trees:  canopy  bottoms 
and  other  dark 
vegetation 

1,955,053  (33.3%) 

4,975,598  (33.1%) 

2 — Dark  Green 

Trees:  inner  canopies 

133,748  (2.3%) 

312,953  (2.1%) 

3 — Dark  Gray 

Road  Asphalt 

561,741  (9.6%) 

1,478,834  (9.8%) 

4 — Light  Gray 

Semi-dark  Concrete  and 
concrete/asphalt  mix 

530,138  (9%) 

1,368,249  (9.1%) 

5 — Middle  Green 

Shrubs  and  many  tree 
tops 

397,829  (6.7%) 

1,035,129  (6.9%) 

6— Gold/Gray 

Darker  Concrete: 
building  and  ground 

216,092  (3.7%) 

495,133  (3.3%) 

7 — Brown 

Dirt,  Mulch,  Sand,  Clay 

655,211  (11.2%) 

1,754,284  (11.7%) 

8 — Bright  Green 

Grass  and  ivy 

442,472  (7.5%) 

1,084,551  (7.2%) 

9— Off  White 

Light  Concrete 

26,916  (0.45%) 

78,452  (0.5%) 

0 — Blue 

Water  and  Bathymetry 
(masked  points) 

950,638  (16%) 

2,439,619  (16.2%) 

Total  Point  Count 

5,869,838 

15,022,802 

*  Capitalized  words  in  the  name  and  description  column  are  the  assigned  class  names. 


1.  K-means  Single  Flight  Line — Central  Campus 

Figure  75  is  a  wide-area  screenshot  of  the  K-means  classifier  results  for  the  single 
flight  line  instance.  It  shows  all  of  the  classes  in  Table  8  other  than  the  water/bathymetry 
class.  Again,  Figure  75  is  the  central  campus  (FL  10  only).  Several  areas  are  labelled  to 
provide  a  sense  of  orientation. 


98 


Figure  75.  K-means  Single  Flight  Line — All  Classes  but  Water. 


From  Figure  75,  it  is  obvious  that  K-means  demonstrates  more  weaknesses  than 
strengths  when  classifying  spectral  LiDAR  data.  Remember,  this  unsupervised  classifier 
does  not  use  the  5%  training  classes  from  N-D  VIS.  Figure  76  shows  the  individual  K- 
means  classes  that  we  loosely  place  into  the  strengths  category  for  the  single  flight  line 
point  cloud.  Quick  Terrain  Modeler’s  “QTA  Discrete  Attribute  Analysis”  tool  allows  one 
to  view  or  hide  points  by  classification,  AGL,  or  any  other  individual  attribute  from 
Table  6. 
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Figure  76.  K-means  Single  Flight  Line  Strengths. 
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The  correct  classification  of  ground  points  along  the  flight  line  edge  represents  the 
primary  benefit  of  running  the  K-means  classifier  on  a  single  flight  line  vice  combined 
flight  lines.  The  majority  of  the  NPS  Quad  in  Figure  76  is  correctly  classified  into  the 
dirt,  mulch,  sand,  and  clay  class.  Additionally,  K-means  does  a  decent  job  distinguishing 
the  interfaces  among  the  road  asphalt,  grass,  and  dirt/mulch/sand/clay  classes:  notice  the 
sharp  definition  of  the  roads  around  Herrmann  Hall  and  the  Arizona  Cactus  Garden’s 
clear  borders  between  dirt/sand  and  grass.  Lastly,  the  light  building  concrete  class  is 
defined  well  by  the  K-means  classifier,  especially  on  the  rooftop  of  King  Hall.  This  is  not 
considered  the  primary  strength  for  single  flight  line  K-means,  because  the  combined 
flight  line  K-means  and  the  ML  classifiers  also  define  the  light  building  concrete  class 
with  exceptional  distinction. 

Nevertheless,  the  K-means  classification  process  generates  serious  problems,  as 
seen  in  Figures  77  and  78.  In  Figure  77,  the  turf  around  the  NPS  flagpole  and  the  adjacent 
asphalt  parking  lot  are  misclassified  as  trees.  Turf  and  areas  of  darker  asphalt  have  low 
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spectral  intensities  for  all  three  wavelengths,  similar  to  tall  and  dense  vegetation  which 
also  experiences  weak  signal  returns.  In  contrast,  tree  tops  are  often  misclassified  as 
shrubs  due  to  the  brighter  returns  from  that  part  of  the  canopy.  Next,  the  rooftops  with 
clay-tile  shingles  (e.g.,  Herrmann  Hall)  demonstrate  gross  misclassifications.  Despite 
having  AGL  as  an  input  metric,  the  K- means  classifier  cannot  separate  dirt,  mulch,  sand, 
and  clay  on  the  ground  from  the  clay  shingles  at  much  greater  heights.  Also,  a  substantial 
amount  of  points  on  the  clay-tile  roofs  are  misclassified  as  semi-dark  concrete  and  some 
are  even  included  in  the  tree  classes.  Furthermore,  K- means  sometimes  has  difficulty 
distinguishing  between  sand  and  concrete,  especially  for  large  areas  of  sand  (e.g., 
Monterey  Beach — not  shown). 


Figure  77.  K-means  Single  Flight  Line  Weaknesses. 


Figure  78  shows  another  K-means  AGL-related  discrepancy,  similar  to  the 
aforementioned  problem  of  ground  dirt,  mulch,  sand,  and  clay  versus  the  clay-tile 
shingles.  The  darker  concrete  (gold)  points  on  the  rooftop  of  Root  Hall  are  grouped  into 
the  same  class  as  some  of  the  adjacent  ground  concrete.  In  Figure  78,  the  trees  (classes  1 


101 


and  2)  and  shrubs  (class  7)  are  hidden,  creating  the  apparent  black  holes  in  the  image. 
Also,  the  entire  rooftop  of  Ingersoll  Hall  is  classified  as  road  asphalt,  however  it  is 
mainly  comprised  of  red  lava  rocks — a  material  class  unidentified  by  K-means. 


Figure  78.  K-means  Single  Flight  Line  Weaknesses  Continued. 


2.  K-means  Combined  Flight  Lines — Entire  Campus 

Running  the  K-means  classifier  on  the  entire  NPS  point  cloud  (four-combined 
flight  lines)  generates  no  significant  improvements  over  the  single  flight  line  instance. 
There  is  slightly  less  confusion  between  sand  and  concrete  in  the  combined  point  cloud. 
Moreover,  the  combined  flight  lines  K-means  demonstrates  all  of  the  serious  weaknesses 
of  the  single  flight  line  instance,  plus  an  additional  issue  in  the  areas  of  flight  line 
overlap.  The  NPS  Quad  in  Figure  79  is  comprised  of  points  from  three  flight  lines  (10, 
12,  and  15),  in  contrast  to  the  single  flight  line  portion  of  the  Quad  in  Figure  76.  Due  to 
variations  in  aircraft  altitude,  speed,  and  other  factors  during  the  flight  campaign,  the 
points  of  different  flight  lines  have  a  slight  vertical  offset  (3-5  cm).  This  becomes  a 
problem  in  the  areas  where  two  or  more  flight  lines  overlap.  K-means  interprets  many  of 
these  offset  points  as  grass,  as  seen  in  the  majority  of  the  Quad  in  Figure  79. 
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Figure  79.  K-means  Combined  Flight  Lines — Additional  Weakness. 


Overall,  K-means  performs  better  than  initially  expected.  Yet,  it  is  insufficient  for 
both  single  flight  line  and  combined  flight  line  classification.  Even  with  the  AGL  metric, 
the  K-means  classifier  fails  to  separate  similar  ground  and  building  materials:  clay-tile  vs. 
dirt,  sand  vs.  concrete,  and  sidewalk  vs.  building  concretes.  The  autonomous  K-means 
classifier  also  struggles  with  spectrally  dark  materials,  including  trees,  turf,  and  asphalt. 
Now,  we  will  demonstrate  how  the  spectral  classification  process  is  improved  by  a 
supervised  classifier  that  operates  based  on  input  classes  from  the  5%  training  subset. 

C.  MAXIMUM  LIKELIHOOD  SUPERVISED  CLASSIFICATION 

ML  is  the  best  classification  method  for  extrapolating  the  N-D  VIS  spectral 
processing  approach,  demonstrated  in  Chapter  IV,  to  the  remaining  unclassified  points. 
First,  the  ML  classifier  creates  a  unique  rule  image  for  each  of  the  selected  input  classes. 
Every  pixel  (laser  data  point)  is  represented  in  each  rule  image  and  assigned  class- 
probability  values  based  on  normally  distributed  class  statistics  (Harris  Geospatial 
Solutions  2017b).  The  ML  algorithm  then  builds  one  classification  image  from  the 
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separate  rule  images.  In  the  classification  image,  every  laser  data  point  is  assigned  a 
classification  that  corresponds  to  the  rule  image  with  the  highest  probability  value  for  that 
point  (2017).  Afterward,  the  rule  images  are  discarded  unless  the  user  chooses  to  save 
them.  According  to  Harris,  ENVI  applies  the  discriminate  ML  formula  described  in 
Figure  80  (Harris  Geospatial  Solutions  2017b). 

Figure  80.  ML  Probability  Equation.  Adapted  from  Richards  (2013,  251). 

Qiix)  =  lnp(cot-)—  V2ln 

g^x)  =  probability  that  data  point  ‘x’  belongs  to  class  ‘i’ 
i  =  class 

x  =  ^-dimensional  data  (where  n  is  the  number  of  bands) 

/?(©,)  =  probability  that  class  co*  occurs  in  the  image  and  is  assumed  the  same  for  all 
classes 

|E,|  =  determinant  of  the  covariance  matrix  of  the  data  in  class  to, 

Er1  =  inverse  of  the  data  covariance  matrix 
an  =  mean  vector 

Figure  81  is  the  probability  graph  for  the  concrete-class  rule  image.  It  somewhat 
resembles  a  fat-tailed  log-normal  distribution,  albeit  with  heavy  skew.  ENVI  applies  a 
logarithmic  scale  to  the  probabilities  on  the  x-axis.  The  probability  of  a  point  being 
concrete  is  zero  on  the  right  edge  of  the  graph  and  increases  toward  one  moving  to  the 
left.  The  green  and  red  lines  in  Figure  81  correspond  to  the  colored  areas  in  Figure  82. 
This  shows  the  points  most  likely  to  be  classified  as  concrete  when  the  classification 
image  is  constructed  from  the  rule  images.  The  points  that  remain  colored  white  in  Figure 
82  are  less  likely  to  be  classified  as  concrete. 
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Number  of  Occurence 


Figure  81.  ML  Probability  Graph  for  Concrete-Class  Rule  Image 


Figure  82.  ML  Points  Colored  According  to  Rule  Graph  for  Concrete  Class. 
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The  ML  probability  threshold  is  an  input  parameter  that  can  be  set  either 
universally  for  all  classes  or  individually  for  each.  While  it  does  not  factor  into  the  ML 
equation,  the  probability  threshold  will  prevent  a  point  from  being  classified  into  the 
most-likely  class  if  that  point  does  not  meet  the  minimum  probability  requirement.  For 
example,  an  arbitrary  point  has  the  following  rule  image  probabilities:  0.1  for  trees,  0.2 
for  asphalt,  0.5  for  grass,  and  0.1  for  turf.  Normally,  this  point  will  get  classified  as  grass 
in  the  classification  image.  However,  if  the  probability  threshold  for  the  grass  class  is  set 
to  0.75,  then  this  point  does  not  meet  the  threshold  and  remains  unclassified — it  is  not 
reviewed  again  for  inclusion  into  any  other  class,  even  the  next  most  likely. 

We  experiment  with  different  probability  thresholds  to  see  if  this  factor  can 
improve  the  results.  We  expect  that  selecting  “no  probability  threshold”  across-the-board 
will  result  in  a  point  cloud  with  many  errors  of  commission  (e.g.,  numerous  points  will  be 
classified  as,  say,  concrete  that  are  not  actually  concrete).  Indeed,  this  is  true.  So,  will 
making  the  classifier  more  discriminate,  by  increasing  the  probability  thresholds,  result  in 
fewer  errors  of  commission  and  more  distinctly  defined  classes?  The  answer:  yes,  there 
will  be  fewer  errors  of  commission,  however  the  cost  is  an  increase  in  errors  of  omission 
(e.g.,  some  of  the  points  that  truly  are  concrete  will  remain  unclassified).  Figure  83  is  a 
comparison  of  ML  results  for  the  clay-tile  shingles  class  at  no  probability  threshold  (left) 
and  50%  probability  threshold  (right). 


Figure  83.  Clay-tile  Shingles  Class  Probability  Thresholds — Errors  of 

Commission  vs.  Omission. 
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Table  9  is  the  point  classification  breakdown  for  the  four  runs  of  the  ML 
classifier:  zero  and  0.5  probability  thresholds  for  both  a  single  flight  line  (FL  10)  and  the 
four-combined  flight  lines  of  the  entire  campus.  Approximate  percentages  are  provided 
for  a  comparison.  A  0.5  probability  threshold  results  in  34%  of  the  point  cloud  remaining 
unclassified — millions  of  points!  Instead  of  attempting  to  find  a  threshold  balance 
between  zero  and  0.5,  we  proceed  with  no  probability  threshold  and  apply  the  post¬ 
processing  technique  in  Section  D.  Post-processing  by  multivariate  filtering  offers  the 
advantage  of  removing  errors  of  commission  without  introducing  errors  of  omission, 
unlike  what  happens  by  simply  increasing  the  probability  threshold  of  the  ML  classifier. 


Table  9.  ML  Point  Breakdown  by  Class  and  Probability  Threshold. 


ML  Single  Flight  Line 

ML  Combined  Flight  Lines 

Class 

No  Threshold 

50%  Threshold 

No  Threshold 

50%  Threshold 

0 — Unclassified 

— 

1,966,674  (34%) 

— 

5,095,287  (34%) 

1 — Grass 

347,778  (6%) 

204,170  (3.4%) 

881,438  (5.9%) 

491,714(5.5%) 

2— Turf 

70,032  (1.2%) 

22,359  (0.4%) 

177,813  (1.2%) 

57,020  (0.38%) 

3 — Clay-tile  Shingles 

104,500  (1.8%) 

31,788  (0.5%) 

254,587(7.7%) 

72,451  (0.48%) 

4 — Dirt  and  Mulch 

999,086  (17%) 

734,721  (12.4%) 

2,594,605 

(17.3%) 

1,880,332 

(72.5%) 

5 — Sidewalk  Concrete 

160,455  (2.7%) 

95,650  (1.6%) 

438,732  (2.9%) 

255,037  (7.7%) 

6 — Red  Brick 

247,467  (4.2%) 

138,904  (2.4%) 

620,056  (4.1%) 

347,399  (2.3%) 

7 — Shrubs 

488,686  (8.3%) 

278,069  (4.6%) 

1,258,132  (8%) 

737,906  (4.9%) 

8 — Trees 

1,263,462 

756,061 

3,178,190 

1,900,236 

(21.5%) 

(12.9%) 

(21.2%) 

(72.6%) 

9 — Dark  Bldg  Concrete 

132,314  (2.3%) 

70,869  (1.1%) 

332,166  (2.2%) 

165,446  (1.1%) 

10 — Tan/Gray  Shingles 

42,053  (0.7%) 

6,042  (0.1%) 

104,734  (0.7%) 

17,830  (0.7%) 

1 1 — Asphalt  Shingles 

81,462(1.4%) 

9,570  (0.2%) 

191,858  (7 .3%) 

17,209  (0.7%) 

12 — Red  Basalt/Lava 
Rock 

60,666  (1%) 

16,725  (0.3%) 

161,088  (7%) 

49,104  (0.3%) 

13 — Sand 

208,664  (3.6%) 

138,864  (2.4%) 

544,199  (3.6%) 

353,910  (2.4%) 

14 — Road  Asphalt 

640,412  (11%) 

404,635  (6.9%) 

1,622,459 

(70.7%) 

1,026,055  (6.8%) 

15 — Light  Bldg 

Concrete 

72,163  (1.2%) 

44,099  (0.8%) 

223,126  (1.5%) 

116,247  (0.8%) 

16 — Water/Bathy 

950,638  (16%) 

950,638  (16%) 

2,439,619 

(76.2%) 

2,439,619 

(76.2%) 

Total  Point  Count 

5,869,838 

15,022,802 
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1.  ML  Single  Flight  Line — Central  Campus 

After  a  few  practice  runs  of  the  ML  classifier,  it  became  apparent  that  some  of  the 
training  classes  in  Table  7  must  be  excluded  from  the  process.  The  general  building  class 
is  too  spectrally  diverse;  it  causes  an  excessive  misclassification  of  vegetation  points  as 
general  building.  Powerlines  have  extremely  low  intensities  in  all  three  wavelengths  (see 
Figure  71).  Consequently,  there  is  excessive  confusion  between  true  powerlines  and  the 
similarly  dark  inner  tree  canopies.  Henceforth,  the  pertinent  ML  input/output  classes  are 
grass  (1),  turf  (2),  clay-tile  shingles  (3),  dirt  and  mulch  (4),  sidewalk  concrete  (5),  red 
brick  (6),  shrubs  (7),  trees  (8),  dark  building  concrete  (9),  tan/gray  shingles  (10),  asphalt 
shingles  (11),  red  basalt/lava  rock  (12),  sand  (13),  road  asphalt  (14),  and  light  building 
concrete  (15).  Again,  water  and  bathymetry  points  (16)  are  carried  over  to  the  classified 
point  cloud  via  masking.  The  following  bands  are  selected  as  input  metrics  for  the  ML 
classifier:  AGL  (Band  11),  Red  (Band  14),  Green  (Band  15),  and  Blue  (Band  16). 
Classification  results  are  not  noticeably  improved  by  including  the  vegetation  indices: 
GNDVI  (Band  18),  GDVI  (Band  19),  and  GRVI  (Band  20). 

Figure  84  shows  the  results  of  the  ML  classifier  for  a  region  near  Herrmann  Hall. 
This  is  the  single  flight  line  instance  (FL  10),  and  the  figure  key  includes  only  the  classes 
that  exist  in  the  chosen  area.  The  advantages  of  ML  over  K-means  are  readily  apparent: 
The  majority  of  rooftop  points  on  Herrmann  Hall  are  correctly  classified  as  clay-tile 
shingles,  not  as  dirt/mulch  or  concrete.  There  exists  minor  confusion  between  clay-tile 
shingles  and  trees — some  points  on  the  clay-tiled  roof  are  misclassified  as  trees,  and 
some  points  in  the  nearby  tree  canopies  are  misclassified  as  clay-tile  shingles. 

Similar  to  K-means,  ML  does  well  in  distinguishing  the  interfaces  among  dirt/m 
ulch,  sidewalk  concrete  and  road  asphalt,  grass,  red  brick,  and  turf.  Nevertheless,  the  red 
brick  and  turf  classes  demonstrate  minor  confusion.  In  Figure  84,  small  sections  of  the 
asphalt  parking  lot  in  front  of  Herrmann  Hall  are  misclassified  as  turf.  Additionally,  it 
appears  that  some  sections  of  sidewalk  concrete  and  the  majority  of  road  paint  (white-, 
red-,  and  yellow-painted  lines)  are  confused  with  red  brick.  We  were  unable  to  separate 
road  paints  from  the  road  asphalt  in  N-D  VIS  (Chapter  IV).  While  it  is  reasonable  to 

assume  that  road  paint  is  spectrally  closer  to  red  brick  than  to  dark  asphalt,  the  process  is 
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complicated  by  the  narrow  margin  of  spectral  separation  between  these  two  classes 
(Refer  to  red  brick  and  asphalt  in  the  ground  materials  spectral  plot — Figure  65). 
Therefore,  we  attribute  this  problem  to  the  analysist’s  inability  to  distinguish  a  road  paint 
training  class,  rather  than  to  the  ML  classifier.  Lastly,  many  of  the  parked  cars  (also 
undefined  in  the  training  subset)  are  misclassified  as  shrubs — the  reason  is  unknown,  but 
cars  are  within  the  same  AGL  range  as  shrubs  (0-17),  far  below  the  mean  AGL  of  trees 
(127). 


Figure  84.  Screenshot  of  Single  Flight  Line  ML  Results  Near  Herrmann  Hall. 
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2.  ML  Combined  Flight  Lines — Entire  Campus 

The  input  classes  and  spectral  bands  (AGL  and  RGB)  for  the  combined  flight 
lines  instance  of  ML  are  identical  to  those  of  the  single  flight  line  instance.  The  strengths 
and  weaknesses  described  previously  for  the  single  flight  line  point  cloud  are  reproduced 
in  the  combined  flight  lines  point  cloud.  Unlike  K-means,  ML  does  not  breakdown  in  the 
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areas  of  flight  line  overlap.  Notice  the  Quad  in  Figure  85;  the  ground  points  are  correctly 
classified  as  dirt/mulch,  not  as  the  grass  seen  in  Figure  79.  ML  keeps  the  points  classified 
as  trees  in  a  single  class,  as  opposed  to  the  three  K-means  classes  that  contain  various 
parts  of  the  canopies  (see  Table  8).  Still,  some  should-be  tree  points  are  misclassified  as 
clay-tile  shingles  or  red  basalt/lava  rock.  Also,  the  beach  sand  directly  touching  the  ocean 
is  often  misclassified  as  red  brick  or  asphalt — wet  sand  is  spectrally  darker  than  dry  sand. 
Similarly,  the  concrete  of  Spanagel  Hall  is  darker  in  color  than  many  of  the  other 
buildings,  exhibiting  a  higher  degree  of  misclassification  as  trees. 


Figure  85.  Raw  ML  Results  for  the  Entire  NPS  Campus — 
Four  Combined  Flight  Lines. 
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Overall,  the  ML  classifier  performs  significantly  better  than  the  K- means 
classifier.  ML  successfully  uses  the  AGL  metric  to  distinguish  between  ground  and 
roofing  concretes  and  between  dirt/mulch  and  clay-tile  shingles.  By  comparison,  the  K- 
means  classifier  seems  to  ignore  the  AGL  input  metric,  frequently  mixing  similar- 
composition  ground  and  building  points  into  the  same  output  classes.  Section  D  will 
apply  post-processing  to  the  point  cloud  in  Figure  85;  compare  this  figure  to  the  refined 
Figure  91.  Look  specifically  at  the  trees,  the  roofs  of  Herrmann  and  Spanagel  Halls,  the 
library,  and  the  NEX. 

Table  10  is  the  confusion  matrix  for  the  entire  campus  point  cloud,  classified  by 
ML  with  no  probability  threshold.  It  compares  the  points  in  the  5%  training  subset  before 
and  after  classification.  Each  class’  greatest  misclassification  (confusion)  is  in  red  text. 
The  bold-type  along  the  diagonal  highlights  individual  class  accuracies — i.e.,  points 
mapped  from  a  particular  training  class  in  Table  7  to  the  correct  ML  output  class  in  Table 
9.  Bold-type  class  names  have  an  accuracy  of  70%  or  better. 

Table  11  provides  a  quick  summary  for  the  confusion  matrix  in  Table  10.  Again, 
these  results  are  only  for  the  ML  combined  flight  lines  instance  with  no  probability 
threshold.  Classes  are  listed  by  descending  order  of  accuracy  in  Table  11.  Potential 
reasons  are  given  for  each  class’  most  frequent  misclassification(s).  Overall  classification 
accuracy  is  75%,  with  10  of  the  14  classes  demonstrating  individual  accuracies  above 
70%.  The  grass  and  turf  classes  outperform  the  others  with  accuracies  of  99.9%  and  95%, 
respectively.  The  worst  performer  is  sidewalk  concrete  at  46%  due  to  a  high  degree  of 
confusion  with  similar  ground  types:  sand,  older  faded  asphalt  roads,  and  roads  paved 
with  a  lighter-color  asphalt/concrete  mixture  (i.e.,  not  fresh  or  dark  black  asphalt).  Lastly, 
approximately  half  of  the  classes  fall  into  the  range  of  moderate  success:  75-89%. 
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Table  10.  Confusion  Matrix — Raw  ML  Results  vs.  the  5%  Training  Subset. 


Grass 

Turf 

Clay-tile 

Dirt/ 

Mulch 

Sidewalk 

Concrete 

Red 

Brick 

Shrubs 

Trees 

Dark 

Bldg 

Concrete 

Tan/gray 

Shingles 

Asphalt 

Shingles 

Red 

Basalt 

Lava 

Rock 

Sand 

Road 

Asphalt 

Light 

Bldg 

Concrete 

Grass 

(26.204) 

26,181 

(99.9%) 

0 

0 

23 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Turf 

(122) 

1 

116 

(95%) 

0 

0 

0 

l 

0 

0 

0 

0 

0 

0 

0 

4 

0 

Clay-tile 

Shingles 

(6,799) 

3 

0 

5,202 

(77%) 

346 

l 

13 

282 

48 

36 

130 

0 

734 

0 

2 

2 

Dirt/Mulch 

(181,865) 

16,293 

4,262 

0 

120,466 

(66%) 

4,417 

17,258 

0 

0 

1 

0 

0 

0 

6,504 

12,664 

0 

Sidewalk 

Concrete 

(31,872) 

56 

14 

0 

1,443 

14,657 

(46%) 

3,770 

0 

0 

0 

0 

0 

0 

7,300 

4,632 

0 

Red  Brick 
(298) 

2 

0 

0 

1 

9 

204 

(68%) 

0 

0 

0 

0 

0 

0 

0 

82 

0 

Shrubs 

(5,995) 

18 

3 

106 

36 

0 

28 

5,108 

(85%) 

301 

50 

58 

178 

94 

0 

7 

8 

Trees 

(190,083) 

0 

0 

2,902 

1 

0 

0 

21,551 

151,666 

(80%) 

1,283 

1,729 

6,978 

3,950 

0 

0 

23 

Dark  Bldg 

Concrete 

(12,078) 

2 

6 

212 

347 

248 

841 

177 

25 

7,521 

(62%) 

453 

10 

27 

56 

150 

2,003 

Tan/gray  Shingles 

(1,313) 

0 

0 

9 

0 

0 

0 

15 

0 

8 

1,098 

(84%) 

181 

2 

0 

0 

0 

Asphalt 

Shingles 

(819) 

0 

0 

0 

0 

0 

0 

16 

3 

1 

32 

767 

(94%) 

0 

0 

0 

0 

Red  Basalt 

Lava  Rock 

(2,588) 

0 

0 

114 

470 

0 

3 

29 

15 

0 

22 

0 

1,935 

(75%) 

0 

0 

0 

Sand 

(14,179) 

0 

0 

0 

228 

496 

716 

0 

0 

0 

0 

0 

0 

12,441 

(88%) 

298 

0 

Road  Asphalt 

(72,356) 

0 

4,066 

0 

526 

915 

5,844 

0 

0 

0 

0 

0 

0 

130 

60,875 

(84%) 

0 

Light  Bldg 

Concrete 

(6,835) 

2 

0 

1 

5 

37 

0 

0 

0 

724 

0 

0 

0 

0 

0 

6,066 

(89%) 

Total  Accuracy 

414,303  of  553,406 
points 

75% 
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Table  11.  Accuracy  Summary  for  Raw  ML  Confusion  Matrix. 


Class  Name  (#) 

Accuracy 
(%)  ' 

Primary 

Misclassification 

Suspected  Reasons  for  Confusion 

Grass  (1) 

99.9 

Dirt/Mulch 

These  two  classes  naturally  exist  in 
extremely  close  proximity.  (Possibly 
some  subpixel  mixing). 

Turf  (2) 

95 

Road  Asphalt 

Dark  in  all  spectral  bands;  both  are 
ground  materials  (AGL  =  0). 

Asphalt 

Shingles  (11) 

94 

Tan/gray  Shingles 

Gray  and  black  asphalt  shingles  are  often 
interspersed  on  the  same  roof. 

Light  Building 
Concrete  (15) 

89 

Dark  Building 
Concrete 

There  are  several  middle  shades  of 
concrete  on  various  roofs  that  arc  not 
separated  in  N-D  VIS. 

Sand  (13) 

88 

Red  Brick. 

Sidewalk  Concrete 

There  are  many  areas  with  a  thin  layer  of 
sand  over  a  concrete  path  (especially  near 
the  beach).  Some  areas  of  wet  sand  are 
spectrally  darker,  resembling  red  brick. 

Shrubs (7) 

85 

Trees 

There  is  greater  confusion  for  vegetation 
points  around  the  cutoff  AGL  between 
shrubs  and  trees.  (AGL  =  ~  1 6). 

Road 

Asphalt  (14) 

84 

Red  Brick.  Turf 

All  arc  spectrally  dark.  Road  paint  left 
unseparated  from  the  asphalt  in  N-D 

VIS — gets  placed  into  red  brick  class. 

Tan/Gray 

Shingles  (10) 

84 

Asphalt  Shingles 

Same  reason  as  asphalt  shingles  (above). 

Trees  (8) 

80 

Shrubs,  Asphalt 
Shingles 

Spectrally  dark  and  many  tree  points  are 
at  similar  AGLs  to  the  asphalt  rooftops. 

Clay-tile 

Shingles  (3) 

77 

Red  Basalt/Lava 

Rock 

Similar  in  color;  no  visible-red  laser  to 
help  distinguish. 

Red  Basalt 

Lava  Rock  (12) 

75 

Dirt/Mulch.  Clay- 
tile  Shingles 

Same  as  previous  reason  for  confusion 
with  clay  shingles.  Confusion  with  dirt/m 
ulch  is  due  to  a  sub-area  of  corrupt  AGLs 
(described  later). 

Red  Brick  (6) 

68 

Road  Asphalt 

No  visible-red  laser;  both  arc  spectrally 
dark.  Almost  all  the  road  paint  within  the 
asphalt  class  is  classified  as  red  brick. 

Dirt/Mulch  (4) 

66 

All  other  ground 
classes,  especially 
Red  Brick 

Large  and  disperse  class  in  N-D  VIS 

Dark  Building 
Concrete  (9) 

62 

Light  Building 
Concrete,  Red  Brick 

Middle  shades  of  concrete  are  not 
distinguishable  in  N-D  VIS,  getting 
lumped  into  the  dark  concrete  class. 

Unsure  on  the  confusion  with  red  brick. 

Sidewalk 

Concrete  (5) 

46 

Sand,  Road  Asphalt 

Concrete  is  made  from  sand,  and  all  three 
classes  have  AGL  =  0  for  all  points. 

Many  of  the  campus  roads  are  more  of  a 
lighter  concrete/asphalt  mixture. 
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D.  POST  PROCESSING:  REFINEMENT  BY  MULTIVARIATE  FILTERING 

Quick  Terrain  Modeler’s  “Multivariate  Filtering”  tool  applies  simple  logical 
operators  to  filter  points.  An  example  is  illustrated  in  Figure  86  for  the  sidewalk  concrete 
class,  filtered  by  number  of  returns.  Most  areas  of  true  sidewalk  concrete  have  only  one 
return,  whereas  the  majority  of  sidewalk  concrete  points  with  more  than  one  return  are 
actually  dirt  or  sand  under  vegetation  (red  points  in  Figure  86).  Of  course,  this  method 
will  incorrectly  filter  some  true-class  points,  however  the  tradeoff  is  a  substantial 
improvement  over  the  raw  ML  results  without  the  introduction  of  a  significant  amount  of 
errors  of  omission.  In  Figure  86,  the  23,250  filtered  points  are  colored  red,  but  this  tool 
also  allows  the  user  to  hide,  cut,  or  reclassify  the  selected  points. 


Figure  86.  Illustration  of  the  Multivariate  Filtering  Process  for 

Sidewalk  Concrete. 
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We  can  apply  multivariate  filtering  for  individual  attributes  or  any  combination  of 
attributes  from  Table  6.  For  this  effort,  we  post-process  only  the  ML  entire  campus  point 
cloud,  limited  to  the  single  most-beneficial  filtering  attribute:  number  of  returns  (#6).  The 
grass  class  is  not  post-processed  due  to  its  raw  classification  accuracy  of  99.9%.  For  the 
other  14  (non-water)  classes,  number-of-retums  filtering  is  applied  in  a  logical  manner 
that  errs  on  the  side  of  caution:  under  correction  is  preferred  to  over  correction.  E.g.,  this 
method  can  decrease  the  confusion  between  single-return  building  materials  and  multi- 
return  vegetation  classes,  however  it  cannot  be  used  to  correct  the  confusion  between 
shrubs  and  trees  (both  multi-return  classes)  or  between  asphalt  shingles  and  dark  building 
concrete  (both  single-return  classes). 

1.  Select  Before-and-After  Figures 

The  upcoming  figures  compare  the  ML  results  before  and  after  post-processing 
for  two  ground  classes  and  all  six  building  classes.  Figures  87  and  88  respectively 
represent  the  before  and  after  for  turf  (2)  and  road  asphalt  (14).  Figures  89  and  90  are  the 
before  and  after  for  clay-tile  shingles  (3),  dark  building  concrete  (9),  tan/gray  shingles 
(10),  asphalt  shingles  (11),  red  basalt/lava  rock  (12),  and  light  building  concrete  (15). 

Figure  87.  Turf  and  Road  Asphalt  Before  Filtering  by  Number  of  Returns. 
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Figure  88.  Turf  and  Road  Asphalt  After  Filtering  by  Number  of  Returns. 


Figure  89.  Building  Materials  Before  Filtering  by  Number  of  Returns. 
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Figure  90.  Building  Materials  After  Filtering  by  Number  of  Returns. 


^  Clay-tile  Shingles 

9  Dark  Bldg  Concrete 

10  Tan/gray  Shingles 

11  Asphalt  Shingles 

12  Red  Basalt/Lava  Rock 


Light  Bldg  Concrete 


2.  NPS  Campus  Refined  by  Number  of  Returns 

In  both  before- and- after  examples,  filtering  by  number  of  returns  greatly 
improves  class  definition  by  removing  the  majority  of  errors  of  commission  without 
introducing  too  many  errors  of  omission.  Figure  91  is  the  refined  entire  campus  point 
cloud.  A  significant  change  can  be  seen  in  the  trees  and  on  the  building  rooftops.  During 
the  multivariate  filtering  process,  two  major  rooftop  misclassifications  are  discovered  and 
corrected:  On  the  NEX,  many  of  the  clay-tile  shingles  and  some  areas  of  the  flat  concrete 
roof  are  misclassified  as  dirt/mulch.  Likewise,  the  red  lava  rock  section  of  the  library  roof 
is  misclassified  as  dirt/mulch.  It  appears  that  these  areas  incorrectly  received  an  AGL 
value  of  zero  during  the  LASheight  script.  Accordingly,  ML  classifies  them  into  the 
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most-appropriate  ground  class  (dirt/mulch) — we  manually  correct  only  these  specific 
errors  (and  absolutely  no  others),  since  the  fault  is  not  attributable  to  the  classifier. 


Figure  91.  Post- Processed  ML  Point  Cloud — Entire  NPS  Campus 

(Four  Combined  Flight  Fines). 


Clay-tile  Shingles  3 

Dirt/Mulch  4 

Sidewalk  Concrete 

Red  Brick  6 

Shrubs  7 

Trees  8 

Dark  Bldg  Concrete _ 9 _ 

Tan/grav  Shinglcdtfll^'51 
V*-  Asphalt  Shingles  1 1 

Red  Basalt/Lava  Rock  12 

SandHHHI 

1  Road  Asphalt  14 

*>  .  Lighi  Bldg  HFH 


Table  12  is  the  confusion  matrix  for  the  post-processed  results  compared  to  the 
original  5%  training  subset.  The  gray-shade  cells  indicate  changes  from  the  raw  MF 
confusion  matrix  in  Table  10.  A  red-text  entry  represents  each  class’  greatest 
misclassification  after  refinement.  Additionally,  a  summary  chart  for  the  post-processing 
confusion  matrix  is  provided  in  Table  13,  listing  each  class’  individual  improvement  and 
the  overall  improvement. 
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Table  12.  Confusion  Matrix — ML  Post-Processed  by  Number  of  Returns  vs.  the  5%  Training  Subset. 


Grass 

Turf 

Clay-tile 

Dirt/Mulch 

Sidewalk 

Concrete 

Red  Brick 

Shrubs 

Trees 

Dark  Bldg 

Concrete 

Tan/gray 

Shingles 

Asphalt 

Shingles 

Red  Basalt 

Lava  Rock 

Sand 

Road 

Asphalt 

Light  Bldg 

Concrete 

Grass 

(26,204) 

26,181 

(99.9%) 

0 

0 

23 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Turf 

(122) 

1 

116 

(95%) 

0 

0 

0 

l 

0 

0 

0 

0 

0 

0 

0 

4 

0 

Clay-tile 

Shingles 

(6,799) 

3 

0 

5,839 

(86%) 

3 

1 

li 

24 

14 

36 

130 

0 

734 

0 

2 

2 

Dirt/Mulch 

(181,865) 

16,293 

615 

0 

146,132 

(80%) 

3,784 

5,786 

0 

0 

1 

0 

0 

0 

6,504 

2,750 

0 

Sidewalk 

Concrete 

(31,872) 

56 

7 

0 

1,450 

14,806 

(46%) 

3,621 

0 

0 

0 

0 

0 

0 

7,300 

4,632 

0 

Red  Brick 
(298) 

2 

0 

0 

1 

9 

204 

(68%) 

0 

0 

0 

0 

0 

0 

0 

82 

0 

Shrubs 

(5,995) 

18 

1 

99 

24 

0 

10 

5,390 

(90%) 

301 

36 

18 

2 

85 

0 

3 

8 

Trees 

(190,083) 

0 

0 

2,229 

1 

0 

0 

21,551 

164,405 

(86%) 

102 

208 

82 

1,487 

0 

0 

18 

Dark  Bldg 

Concrete 

(12,078) 

2 

4 

212 

286 

244 

810 

66 

21 

7,934 

(66%) 

453 

10 

27 

0 

6 

2,003 

Tan/gray  Shingles 
(1,313) 

0 

0 

9 

0 

0 

0 

3 

0 

8 

1,110 

(85%) 

181 

2 

0 

0 

0 

Asphalt 

Shingles 

(819) 

0 

0 

0 

0 

0 

0 

11 

1 

1 

32 

774 

(95%) 

0 

0 

0 

0 

Red  Basalt 

Lava  Rock 
(2,588) 

0 

0 

114 

3 

0 

2 

9 

2 

0 

22 

0 

2,436 

(94%) 

0 

0 

0 

Sand 

(14,179) 

0 

0 

0 

228 

496 

0 

0 

0 

0 

0 

0 

0 

13,455 

(95%) 

0 

0 

Road  Asphalt 
(72,356) 

0 

3,719 

0 

873 

915 

5,844 

0 

0 

0 

0 

0 

0 

130 

60,875 

(84%) 

0 

Light  Bldg 

Concrete 

(6,835) 

2 

0 

1 

5 

37 

0 

0 

0 

724 

0 

0 

0 

0 

0 

6,066 

(89%) 

Total  Accuracy 

455,723  of  553,406 
points 

82% 
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Table  13.  Accuracy  Summary  for  Post-Processing  ML  Confusion  Matrix. 


Class  Name  (#) 

Raw  ML  Accuracy  (%) 

from  Table  1 1 

Accuracy  After  Post- 
Processing  (%) 

Rate  of 

Improvement  (%) 

Grass (1) 

99.9 

N/A 

N/A 

Turf  (2) 

95 

95 

- 

Asphalt 

Shingles  (11) 

94 

95 

+  1 

Light  Building 
Concrete  (15) 

89 

89 

- 

Sand  (13) 

88 

95 

+7 

Shrubs  (7) 

85 

90 

+5 

Road 

Asphalt  (14) 

84 

84 

- 

Tan/Gray 
Shingles  (10) 

84 

85 

+  1 

Trees (8) 

80 

86 

+6 

Clay-tile 
Shingles  (3) 

77 

86 

+9 

Red  Basalt 

Lava  Rock 
(12) 

75 

94 

+19 

Red  Brick  (6) 

68 

68 

- 

Dirt/Mulch  (4) 

66 

80 

+14 

Dark  Building 
Concrete  (9) 

62 

66 

+4 

Sidewalk 
Concrete  (5) 

46 

46 

- 

Overall 

75 

82 

+7 

Considerable  individual  class  improvements  are  in  bold  type. 


Post-processing  by  number  of  returns  increases  overall  classification  accuracy  to 
82%.  The  dirt/mulch  and  red  basalt/lava  rock  classes  demonstrate  the  greatest  progress  of 
+14%  and  +19%,  respectively.  Areas  of  dirt  and  mulch  now  exhibit  better  distinction 
from  turf  and  the  impervious  surfaces,  primarily  road  asphalt  and  red  brick.  The 
improvement  in  the  red  basalt/lava  rock  class  is  mostly  attributable  to  the  aforementioned 
correction  of  the  library  roof — misclassified  as  dirt/mulch  due  to  false  AGL  values.  Had 
this  error  not  occurred,  we  would  expect  red  basalt/lava  rocks  to  demonstrate 
approximately  90%  or  better  classification  accuracy.  Some  confusion  persists  among 
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clay-tile  and  red  basalt/lava  rock  roofs,  which  is  uncorrectable  without  analyzing  other 
attributes  besides  number  of  returns  (or  outright  manual  correction). 

Figure  92  is  the  area  near  Herrmann  Hall  previously  displayed  in  Figure  84,  but 
now  post-processed  by  number  of  returns.  Notice  the  sharper  class  definitions,  especially 
on  the  rooftop;  there  remains  no  noticeable  confusion  among  the  building  materials  and 
vegetation  classes.  The  misclassifications  of  road  paint  as  red  brick,  asphalt  as  turf,  and 
cars  as  shrubs  remain  the  three  apparent  discrepancies  in  the  scene,  however  these  are 
considered  minor — and  uncorrectable  without  additional  post-processing.  Figure  92, 
refined  spectral  LiDAR  classification,  should  also  be  compared  to  Figure  49,  the  best  of 
spatial-only  ( traditional )  classification  plus  manual  correction.  The  increased  accuracy 
and  diversity  in  terrain,  building,  and  vegetation  classification  due  to  the  spectral 
component  is  evident:  82%  (minimally  refined  spectral  accuracy)  vs.  77.36%  (heavily 
refined  spatial  accuracy — from  Table  4)  and  16  classes  vs.  four. 


Figure  92.  Post-Processed  ML  Results  Near  Herrmann  Hall. 
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After  multivariate  filtering,  the  tree  and  shrub  classes  are  unconfused  with  the 
building  materials,  as  seen  in  Figures  91  and  92.  Less  noticeable  is  the  7%  improvement 
for  sand,  primarily  from  the  de-confliction  of  wet  beach  sand  with  asphalt  and  red  brick. 
Eight  of  the  14  classes  maintain  essentially  the  same  self-mapping  accuracies:  grass  (not 
post-processed),  turf,  asphalt  shingles,  light  building  concrete,  road  asphalt,  tan/gray 
shingles,  red  brick,  and  sidewalk  concrete  (still  the  worst  at  46%).  For  these  classes,  any 
improvement  in  class  distinction  comes  primarily  from  the  de-confliction  of  true-class 
points  with  the  dirt/mulch  and  vegetation  classes,  as  exemplified  in  the  before-and-after 
figures  for  turf  and  road  asphalt  (Figures  87  and  88).  Altogether,  post-processing  by 
number  of  returns  makes  a  worthwhile  contribution  to  the  spectral  FiDAR  terrain 
classification  process.  Additional  filtering  attributes  that  may  further  refine  the  results  are 
the  vegetation  indices  (GNDVI,  GDVI,  and  GRVI),  AGF,  and  possibly  RGB.  We  leave 
this  as  future  work  for  the  next  analyst. 
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VI.  CONCLUSION 


A.  APPLICABILITY  TO  SPACE  AND  MILITARY  OPERATIONS 

LiDAR  technologies  continue  to  climb  increasingly  upward  in  operational 
altitude.  Eventually,  a  narrow-beam  multi-wavelength  LiDAR  system  will  become  an 
important  remote- sensing  satellite  payload.  The  analysis  techniques  applied  herein 
represent  a  demonstration-of-concept  that  multi-wavelength  LiDAR  systems  improve 
upon  terrain  classification  by  adding  a  spectral  component.  Lurthermore,  ENVI’s  n- 
Dimension  Visualizer  has  successfully  exhibited  the  capacity  to  adapt  to  LiDAR  data. 
This  powerful  tool  can,  in  theory,  support  spectral  LiDAR  data  of  up  to  “n”  wavelengths 
(i.e.,  a  hyperspectral  LiDAR — currently  an  infant  laboratory  technology).  With  a  greater 
number  of  available  dimensions  (wavelengths),  class  clusters  will  become  even  more 
recognizable  and  easier  to  distinguish  in  spectral  space. 

As  demonstrated,  spectral  LiDAR  imaging  allows  us  to  analyze  3D  spatial  and 
spectral  attributes  simultaneously.  This  may  provide  significant  advantages  for 
identifying  adversarial  deception  tactics,  such  as  camouflage.  By  analyzing  both  the 
spatial  dimensions  and  the  spectral  returns  of  a  scene,  camouflage  is  more  likely  to 
become  distinguishable  from  real  vegetation.  Analysis  techniques  for  multispectral  and 
hyperspectral  electro-optical  (EO)  imagery  also  possess  the  capacity  to  separate  fake 
from  actual  vegetation,  however  LiDAR  is  an  active  remote-sensing  method  that  can  be 
employed  at  night.  Moreover,  it  offers  high-resolution  3D  data  products  that  compliment 
traditional,  lower-resolution  and  2D,  multi-  and  hyperspectral  imagery. 

In  the  military  space  domain,  technological  innovation,  evolution,  and 
adaptability  permit  us  to  remain  ahead  of  our  advancing  adversaries.  Considerable  effort 
should  remain  devoted  to  the  evolution  of  spectral  LiDAR  systems  and  the  maturation  of 
accompanying  data-analysis  techniques.  This  technology  may  prove  tremendously 
beneficial  if  it  becomes  available  from  space  over  areas  denied  even  to  airborne 
intelligence,  surveillance,  and  reconnaissance  (ISR)  platforms. 
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B.  FUTURE  WORK  IN  SPECTRAU  UIDAR 

One  potential  avenue  of  future  work  is  to  analyze  the  spectral  RGB  attributes  in 
Cloud  Compare  as  XYZ  to  identify  class  clusters.  Once  the  flight  line  LAS  files  are 
converted  into  ASCII  files  via  the  nearest  neighbor  script,  the  RGB  attributes  can  be 
imported  into  Cloud  Compare  software  as  x,  y,  and  z  values.  This  creates  a  3D  scatterplot 
similar  to  viewing  the  RGB  in  N-D  VIS.  Spectral  class  clusters  are  selected  and  exported 
as  new  data  classifications.  A  future  effort  can  conduct  the  spectral  analysis  using  Cloud 
Compare  and  then  compare  the  output  to  the  results  from  N-D  VIS  for  the  same  scene. 

Additionally,  an  effort  should  be  made  to  consolidate  the  data  preparations 
process,  discussed  in  Chapter  III,  into  a  single  script.  This  script  can  streamline  the 
process  by  taking  in  the  full-flight  line  LAS  files  (as-received  from  NCALM)  and 
outputting  the  data  in  ENVI  file  format  ready  for  immediate  analysis  in  N-D  VIS. 
Alternately,  the  consolidated  script  can  generate  all  of  the  intermediate  outputs,  so  that 
the  analyst  may  still  review  the  direct  outputs  of  LASground,  LASheight,  LASclassify, 
and  nearest  neighbor. 

Finally,  the  analyst  should  setup  various  camouflaged  targets  in  the  scene  of  a 
future  Optech  Titan  flight  campaign.  Ideally,  military-grade  camouflage  should  be  used  if 
available.  The  N-D  VIS  technique  can  then  be  applied  to  spectrally  separate  the 
camouflage  from  the  true  vegetation.  This  effort  would  have  the  potential  to  demonstrate 
the  advantages  of  spectral-LiDAR  classification  over  traditional  geometric-based  (purely 
spatial)  classifiers  in  a  military-relevant  scenario. 
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APPENDIX.  FOUNDATIONAL  DATA  PROCESSING:  NASA 
DATASET  DERIVED  GRAPHS,  POINT  CLOUDS,  AND  DIGITAL 

ELEVATION  MODELS 


The  purposes  of  this  appendix  are  to  outline  preliminary  work  with  laser  point 
cloud  modeling  and  processing  software  (Sections  B  and  C)  and  to  illustrate  that  multi¬ 
wavelength  LiDAR,  whether  for  topographic  or  atmospheric  purposes,  provides  more 
information  than  monochromatic  LiDAR  alone  (Section  A).  This  represents  important 
groundwork  for  the  spectral  LiDAR  classification  techniques  demonstrated  in  Chapters 
TTT-V.  Quick  Terrain  Modeler  (QTM)  and  Cloud  Compare  (CC)  are  utilized  to  develop 
conventional  LiDAR  datasets  (in  ASCII  and  LAS  formats)  from  two  archived  NASA 
laser  altimetry  data  records:  the  Mars  Orbiter  Laser  Altimeter  (MOL A)  and  the 
Clementine  Lunar  Probe.  Additionally,  atmospheric  graphs  are  produced  for  a  single 
mission  day  of  the  Phoenix  Mars  Lander’s  meteorological  LiDAR  station. 

The  MOLA  and  Clementine  datasets  are  converted  into  spreadsheet  format  for 
preparatory  work  in  Microsoft  Excel,  then  imported  into  QTM  and  CC  to  create  two- 
dimensional  (2D)  and  three-dimensional  (3D)  point  clouds  of  Mars  and  the  Moon.  For 
the  MOLA  dataset,  some  additional  orbital  graphs  are  made  in  Excel.  The  exact 
processes  are  described  thoroughly  in  Section  B  for  MOLA  and  Section  C  for 
Clementine.  Figure  93  is  a  general  workflow  diagram  for  both  datasets. 

The  Phoenix  work  in  Section  A  does  not  follow  Figure  93,  as  it  is  constructed 
using  Excel  alone.  The  Phoenix  atmospheric  graphs  reveal  how  the  combination  of  532 
nm  and  1064  nm  laser  information  generates  a  greater  amount  of  understanding  about  the 
composition  of  the  Martian  atmosphere.  Although  the  bulk  of  this  appendix  involves 
single-,  not  multi-wavelength,  LiDAR,  it  stands  as  a  necessary  learning  experience  and 
foundation  for  the  main  thesis  effort. 


129 


Figure  93.  Workflow  Diagram  for  NASA  Laser  Altimetry  Datasets. 


A.  PHOENIX  MARS  LANDER  “PHX-M-MET-3-L-RDR- V  1.0” 

The  NASA  higher-level  Reduced  Data  Record  (RDR)  is  selected  over  the  raw 
dataset,  the  Experiment  Data  Record  (EDR),  to  create  graphs  from  Phoenix  atmospheric 
LiDAR  data.  The  RDR  has  a  more  intuitive  data  format  than  the  EDR.  For  example,  the 
RDR  gives  532  nm  and  1064  nm  analog  return  signals  in  volts,  whereas  the  EDR’s 
LiDAR  return  signals  are  in  the  form  of  digital  numbers  that  would  require  a  conversion 
to  volts.  This  dataset  is  available  via  file  transfer  protocol  (FTP)  download  from  the 
Washington  University  (St.  Louis,  Missouri)  Node  of  the  NASA  Planetary  Data  System 
(PDS ) ,  http ://pds -geoscience s .  wu stl . edu/default .htm. 

1.  Martian  Atmospheric  Graphs  for  Sol  100  (05  Sep  2008) 

For  this  effort,  a  one  day  subset  of  Phoenix  Lander  data  is  selected  to  explain  the 
concept  that  multi-wavelength  LiDAR  provides  additional  useful  information,  even  for 
places  without  vegetation,  buildings,  or  large  lifeforms  (i.e.,  Mars).  The  lander’s 
meteorological  station  collected  the  LiDAR  data  for  Sol  100  over  a  period  of 

approximately  one  hour,  from  23:52  on  5  September  2008  to  00:52  on  6  September  2008. 
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The  aforementioned  are  the  local  mean  solar  times  at  the  lander’s  location  near  the 
Martian  north  pole:  68.218°  north  latitude,  234.248°  east  longitude. 

The  following  graphs  are  Martian  atmospheric  profiles  as  a  function  of  return 
signal  voltage  for  two  distinct  moments  in  time:  the  first  and  last  cycles  of  the  Sol  100 
measurement  period.  For  laser  backscatter  measurements  below  10  km,  Phoenix  uses 
analog  detection  for  both  wavelengths:  532  nm  and  1064  nm  (Dickinson  et  al.  2008). 
Photon  counting  is  employed  for  532  nm  backscatter  detection  from  between  10  km  and 
20  km;  no  method  exists  on  Phoenix  for  1064  nm  backscatter  detection  from  altitudes 
above  10  km  (Dickinson  et  al.  2008).  Here,  only  the  analog  signals  for  backscatter  returns 
below  10  km  are  graphed.  Figure  94  corresponds  to  the  1064  nm  (infrared)  return  signal, 
and  Figure  95  corresponds  to  the  532  nm  (green)  return  signal. 


Figure  94.  Phoenix  Mars  Lander,  Sol  100:  1064  nm  LiDAR  Backscatter  vs. 

Altitude.  Adapted  from  Dickinson  et  al.  (2008). 
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In  both  figures,  the  y-axes  display  a  logarithmic  scale  of  the  altitude  above  the 
Phoenix  Lander,  from  0  m  to  10,000  m.  The  x-axes  are  ranges  for  average  return  signal 
voltage  (due  to  backscatter  from  the  atmosphere).  LiDAR  backscatter  voltage  is  recorded 
in  vertical  bins  of  20  m,  so  a  data  point  exists  for  every  20  m  of  altitude.  Due  to  the 
logarithmic  scale  of  the  y-axes,  the  data  points  above  1,000  m  converge  to  resemble  a 
solid  line. 


Figure  95.  Phoenix  Mars  Lander,  Sol  100:  532  nm  LiDAR  Backscatter  vs. 

Altitude.  Adapted  from  Dickinson  et  al.  (2008). 
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The  depicted  Sol  100  profiles  suggest  concentrations  of  atmospheric  dust  between 
the  altitudes  of  100  m  and  800  m,  with  the  overall  volume  of  dust  increasing  from  23:52 
to  00:52.  On  this  night,  it  seems  that  no  Martian  clouds  exist  above  the  dust  layer,  as 
indicated  by  the  straight  vertical  lines  from  approximately  1  km  to  10  km.  This 
conclusion  of  an  increase  in  atmospheric  dust  over  the  hour  is  further  supported  by  the 
greater  backscatter  of  the  532  nm  signal  at  00:52,  as  compared  to  the  1064  nm  return 
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signal  at  the  same  time.  The  shorter  wavelength  makes  the  532  nm  photons  more 
susceptible  to  backscattering  from  atmospheric  dust  particles.  Here,  multi-wavelength 
LiDAR  offers  us  more  information  about  the  Martian  atmosphere  than  monochromatic 
LiDAR  alone.  In  “Clouds  and  Precipitation  at  the  Phoenix  Mars  Lander  Site,”  Dickson  et 
al.  describe  Sol  100  as  “dust  only:  exhibiting  a  smooth  profile  of  extinction  coefficients” 
(Dickson  et  al.  2012,  1). 

B.  MARS  ORBITER  LASER  ALTIMETER  “PRECISION  EXPERIMENT 

DATA  RECORD” 

The  NASA  Precision  Experiment  Data  Records  (PEDRs)  are  selected  over  the 
Aggregated  Experiment  Data  Records  (AEDR)  to  save  time  in  the  production  of  high 
quality  2D  and  3D  terrestrial  LiDAR  products  of  Mars.  The  PEDR  dataset,  officially 
identified  as  “MGS-M-MOLA-3-PEDR-L1A-V1.0,”  is  the  result  of  precision  orbital, 
geometric,  and  equipment  calibration  corrections  applied  by  NASA  to  the  raw  AEDR 
files  (Smith  et  al.  1999).  Released  in  1999,  the  dataset  is  available  via  FTP  from  the 
Washington  University  (St.  Louis,  Missouri)  Geosciences  Node  of  the  NASA  PDS, 
http://pds-geosciences.wustl.edu/default.htm. 

Initially,  the  771  PEDR  files  are  downloaded  as  compressed  binary  files  (e.g., 
“ap015781.b”).  The  number  within  each  filename  roughly  corresponds  to  the  beginning 
orbit  of  the  file.  Each  PEDR  binary  file  contains  data  for  several  MOLA  orbits, 
approximately  500,000  to  750,000  laser  shots  per  file.  The  NASA-programmed 
command-line  script  “PEDR2TAB”  is  used  to  extract  the  data  as  TAB  files,  viewable  in 
Microsoft  Excel.  This  program  is  courtesy  of  MOLA  Team’s  Gregory  Neumann,  and  it 
includes  a  version  for  Windows,  MAC  OS,  and  Sun.  The  accompanying  FORTRAN 
source  code  file  “PEDR2TAB.PRM”  contains  several  lines  of  data  extraction  options.  A 
text  editor  places  a  “T”  on  each  line  desired  for  output  to  the  TAB  file.  For  this  effort, 
PEDR2TAB.PRM  line  option  “1”  and  “ground  returns  only”  are  selected  to  extract  the 
following  attributes  for  each  laser  point:  IAU  2000  longitude  (east  positive),  IAU  2000 
latitude  (north  positive),  topography  (planetary  radius  minus  areoid),  MOLA  range, 
planetary  radius,  trigger  channel,  attitude  flag,  shot  number,  packet  sequence  number,  and 
orbit  number.  Figure  96,  shows  a  portion  of  an  extracted  PEDR  file. 
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Figure  96.  Sample  from  Extracted  MOLA  PEDR  File. 

Source:  Smith  et  al.  (1999). 
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1.  Processing  of  MOLA  PEDR  Files  to  Create  Data  Products 

To  improve  compatibility  with  Cloud  Compare  “3D  point  cloud  and  mesh 
processing  software,”  all  of  the  771  MOLA  TAB  files  are  opened  individually  in 
Microsoft  Excel  and  converted  into  comma  delimited  text  format  (.csv).  Additionally,  the 
non-essential  data  columns  are  removed  from  each  file  (see  Figure  96)  in  order  to  reduce 
the  size  of  the  final  dataset. 

After  initial  processing,  each  .csv  file  contains  four  columns  of  data:  longitude, 
latitude,  planetary  radius,  and  topography.  Longitudes  and  latitudes  are  in  IAU  2000 
format;  planetary  radii  and  topographies  are  in  meters.  Planetary  radii  represent  the 
distances  from  the  NASA-determined  gravitational  center  of  Mars  to  the  location  of  the 
MOLA  laser  hit  on  the  surface.  The  dataset  documentation  explains  that  NASA 
calculates  planetary  radii  from  an  areoid  of  stable  rotational  and  gravitational  potential, 
based  on  the  Goddard  Mars  Gravity  Model:  mgm  1025  (Smith  et  al.,  1999  and  2003). 
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2.  2D  Process  and  Results — IAU  2000  (Lat/Long)  Coordinates 
a.  MOLA  Topography  Graphs  for  Orbit  #  20328 

First,  two  graphs  are  produced  from  the  .csv  file  containing  data  for  orbit  number 
20,328.  Figure  97  is  a  single-orbit  representation  of  the  Mars  Global  Surveyors’  ground 
track  derived  from  actual  MOLA  ground  returns.  The  x-axis  is  in  IAU  2000  longitude  (0- 
360°  East),  and  the  y-axis  is  in  IAU  2000  latitude  (-90-90°  North).  The  large  gap  at 
approximately  170°  East  most  likely  results  from  the  decision  to  extract  only  valid  ground 
returns  via  the  aforementioned  PEDR2TAB  option:  “1 — Ground  Returns  Only.’ 

Additionally,  the  longitudes  and  the  corresponding  topography  values  are 
combined  in  Figure  98  to  make  an  Apollo-style  topography  graph  (See  Figure  2)  for 
Mars.  Several  prominent  Martian  landscape  features  are  identified  in  red  with  the 
longitude  and  latitude  coordinates  in  parenthesis.  Again,  the  vertical  line  at  170°  East 
corresponds  to  the  gap  in  valid  ground  hits  also  seen  in  Figure  97. 


Figure  97.  Orbit  Track  Based  on  MOFA  Ground  Hits — Orbit  20,328.  Adapted 

from  Smith  et  al.  (1999). 
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Figure  98.  MOLA  Topography  Above  Mean  Radius — Orbit  20,328.  Adapted 

from  Smith  et  al.  (1999). 


b.  2D  Maps  and  DEM 

All  of  the  771  .csv  files  are  loaded  in  succession  into  Cloud  Compare  v2.6.0  using 
longitudes  as  the  x-values,  latitudes  as  the  y-values,  and  topographies  as  the  z-values. 
Next,  the  minimum  and  maximum  z-values  for  the  771  Cloud  Compare  entities  are 
inspected  to  remove  values  outside  of  valid  Martian  topographies:  -9,200  m  to  21,274  m. 
“Apl8004.csv”  is  discarded  entirely,  because  it  contains  a  large  number  of  invalid 
topography  values.  A  few  additional  .csv  files  contained  one  or  two  invalid  topographies; 
we  found  and  removed  these  errors  quickly  via  basic  Excel  sort  and  filter  functions. 
Seven  hundred  seventy  files  are  then  reloaded  into  Cloud  Compare.  These  entities  are 
now  ready  for  merging. 

The  Cloud  Compare  “Merge  Multiple  Clouds”  function  compiles  the  770  separate 
point  clouds  into  one  entity  of  approximately  590  million  points.  We  export  this  to  ASCII 
immediately.  The  resulting  combined  point  cloud,  “MOLA  2D  ASCII.txt,”  occupies  29.7 
gigabytes  of  disk  space.  Figure  99  is  a  screenshot  of  the  MOLA  2D  point  cloud  opened  in 
Quick  Terrain  Modeler  v8.0.5.2  with  elevations  above  areoid  displayed  in  grayscale. 
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Figure  99.  Ungndded  MOLA  2D  Point  Cloud  (QTM) — 
Topography  in  Grayscale. 


A  generic  geodetic  coordinate  system  is  necessary  in  QTM  to  view  the  point 
cloud  properly  in  the  3D  viewer  mode,  because  degrees  of  longitude  and  latitude  are  the 
x-  and  y-values,  respectively.  The  default  Cartesian  coordinate  system  projects  the  map 
(Figure  99)  properly  in  the  2D  viewer  mode,  but  not  in  the  3D  mode. 

To  produce  the  complete  planetary  DEM  in  Figure  100,  the  point  cloud  is 
converted  to  a  QTT  (Quick  Terrain  Modeler  Gridded  Surface)  using  the  “Convert  Model 
Type  Function.”  We  accept  the  following  default  options:  grid  sampling  at  -0.01, 
adaptive  triangulation  fill,  and  max  z  algorithm.  Larger  grid-sampling  settings  result  in  a 
gridded  surface  that  distorts  the  planet’s  natural  features,  and  settings  less  than  0.01  result 
in  program  crashes  due  to  the  large  data  size  of  the  point  cloud.  Based  on  the  good  result 
from  the  default  options  (Figure  100)  and  our  disappointing  tests  with  different  grid 
sampling  values,  we  decided  not  to  experiment  with  altering  the  other  parameters. 


137 


Figure  100.  MOLA  Digital  Elevation  Model  of  Mars  (QTM) — 

Topography  in  Color. 


Figures  99  and  100  maintain  the  IAU  2000  coordinate  system  and  display  it  in  a 
cylindrical  projection.  Zero  degrees  longitude  corresponds  to  the  left  map  edge,  and  359 
degrees  longitude  corresponds  to  the  right  edge.  NASA-produced  DEMs,  included  in  the 
Mission  Experiment  Gridded  Data  Records  (MEGR),  have  resolutions  of  approximately 
128  pixels  per  degree  of  longitude  at  the  equator  and  mid-latitudes  and  512  pixels  per 
degree  at  the  poles — pixel  size  is  approximately  0.463  km/pixel  (Smith  et  al.  2003).  The 
resolutions  of  the  figures  in  this  section  exhibit  negligible  degradation  from  the  NASA 
products,  even  with  the  one  discarded  file,  “apl8004.csv.”  The  data  loss  is  approximately 
700,000  of  590,000,000  points  (0.0012%). 

The  complete-planetary  ASCII  file  loads  slowly,  and  it  is  difficult  to  navigate  and 
analyze  on  a  computer  with  limited  RAM  (<12  gigabytes).  Thirty  two  gigabytes  or 
greater  RAM  should  alleviate  this  problem.  Quick  Terrain  Modeler’s  “rectangular 
selection”  tool  allows  one  to  select  and  export  a  subset  of  data  points  in  laser  file  format 
(.las).  This  permits  closer  terrain  analysis  without  lengthy  screen  refresh  delays  due  to  the 
large  data  file.  The  area  selected  for  “MOLA  Portion.las”  includes  prominent  Martian 
terrain  features  in  the  western  hemisphere,  such  as  Olympus  Mons,  the  other  Tharsis 
volcanoes,  and  Valles  Marineris. 
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Figure  101  is  the  un-gridded  point  cloud  for  the  “MOLA  Portion”  subset.  It  shows 
the  Mars  Global  Surveyor’s  numerous  orbital  ground  tracks  across  the  landscape.  Figure 
102  is  the  same  scene  in  gridded  format.  Again,  the  default  QTM  “Convert  Model  Type” 
options  are  selected  to  generate  the  rasterized  surface  model.  Remember,  a  geodetic 
coordinate  system  is  required  to  properly  view  the  point  cloud  in  the  3D  viewer  mode  of 
QTM  with  IAU  2000  coordinates  as  x-  and  y-values.  QTM  offers  a  coordinate  conversion 
utility  that  is  straightforward. 


Figure  101.  MOLA  Subset  Point  Cloud  (QTM) — Topography  in  Color. 
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Figure  102.  MOLA  Subset  DEM  (QTM) — Topography  in  Color. 


3.  3D  Conversion  and  Results — IAU  2000  to  Cartesian  Coordinates 

Unfortunately,  no  Cloud  Compare  or  Quick  Terrain  Modeler  functions  exist  that 
can  directly  convert  the  complete  2D  point  cloud  (Figure  99)  into  a  spherical 
representation  of  Mars.  Nevertheless,  it  is  possible  to  derive  the  3D  Cartesian  coordinates 
for  each  laser  point  from  the  NASA-provided  spherical  coordinates:  longitude,  latitude, 
and  planetary  radius.  Table  14  lists  the  parameters  and  formulas  for  the  spherical-to- 
Cartesian  conversion,  applied  individually  for  each  MOLA  point  in  all  770  .csv  files. 
MACRO  functions  in  Excel  are  used  to  rapidly  apply  the  three  formulas  to  each  line  item 
(laser  point),  on  average  500,000  to  750,000  lines  per  file.  To  verify  that  each  file  is 
processed  correctly  via  the  MACROs,  no  batch  processing  is  implemented:  all  770  files 
are  opened,  processed,  inspected,  and  saved  individually. 
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Table  14.  Generic  Mathematical  Formulas  for  Converting  Each  MOLA  Point. 


Spherical  Coordinates 

Polar  Angle  (0) 

=  (90  degrees  -  IAU2000  Latitude) 

Azimuth  Angle  (cp) 

=  IAU2000  Longitude 

r  (meters) 

=  Planetary  Radius 

Derived  Cartesians 

X  (meters) 

=  cos(cp)  x  sin(9)  x  r 

Y  (meters) 

=  sin(cp)  x  sin(9)  x  r 

Z  (meters) 

=  cos(0) x  r 

Now,  each  MOLA  .csv  file  contains  seven  columns:  IAU  2000  longitude,  IAU 
2000  latitude,  planetary  radius,  Cartesian  x,  Cartesian  y,  Cartesian  z,  and  topography.  The 
770  files  are  imported  into  Cloud  Compare  using  the  Cartesian  coordinates  with 
topography  as  a  scalar  value.  The  separate  point  clouds  are  merged  into  one  entity,  and 
the  result  is  saved  as  “MOLA  3D.csv”  (36.2  gigabytes). 

Lor  viewing  the  3D  point  cloud  in  QTM,  it  is  essential  to  import  the  topographies 
as  the  “alpha”  values  and  to  set  the  height  scale  to  “color  by  alpha.”  Also,  the  3D  point 
cloud  must  be  viewed  in  a  Cartesian  coordinate  system,  unlike  the  2D  cloud  that  requires 
a  geodetic  system.  Lastly,  QTM  gridding  cannot  successfully  produce  a  3D  (spherical) 
raster  of  Mars — the  attempt  resulted  in  a  complete  failure  which  resembled  an  accordion 
ball.  However,  the  spatial  point  density  of  the  un-gridded  cloud  is  sufficient  to  project  the 
appearance  of  a  complete  surface,  unless  the  viewer  examines  extremely  close  to  the 
planet — where  you  will  see  ground  tracks  similar  to  Ligure  101. 

Again,  resolution  is  approximately  128  pixels/degree  of  longitude  at  the  equator 
and  512  pixels/degree  at  the  poles.  When  the  display  screen  is  zoomed  in  close  to  the 
planet’s  surface,  the  satellite’s  ground  tracks  (represented  by  crisscrossing  lines  of  laser 
points)  are  significantly  closer  together  at  the  poles  than  at  the  equator.  Ligures  103,  104, 
105,  and  106  are  select  screenshots  of  the  3D  MOLA  point  cloud  in  Cloud  Compare  and 
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Quick  Terrain  Modeler.  The  initials  of  the  respective  program  are  indicated  in 
parentheses  in  the  figure  titles.  Additionally,  each  image  is  color-scaled  differently  to  best 
highlight  the  topographic  features  of  the  particular  region  in  view. 


Figure  103.  MOLA  3D  North  Pole — Topography  in  Color  (CC). 


142 


Figure  104.  MOLA  3D  South  Pole — Color  (QTM). 


Figure  105.  MOLA  3D  Eastern  Hemisphere — Grayscale  (QTM). 
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Figure  106.  MOLA  3D  Western  Hemisphere 
Topography  in  Grayscale  (QTM). 


C.  CLEMENTINE  LUNAR  PROBE  “CLEM1-L-LIDAR-3-TOPO-VLO” 

To  produce  a  Clementine -data  point  cloud  of  the  Moon,  the  “CLEM1-L-LIDAR- 
3-TOPO-Vl.O”  dataset  is  downloaded  via  FTP  from  the  PDS  Geosciences  Node, 
Washington  University  (St.  Louis,  Missouri).  The  dataset  includes  four  TAB  files  plus 
the  associated  label  files,  each  containing  data  for  approximately  100  orbits: 
r009_099.tab,  rl00_199.tab,  r200_299.tab,  and  r300_346.tab.  These  files  are  combined 
into  one  .csv  file  of  over  100,000  laser  data  points.  Next,  each  line  item  is  compared  to 
the  lunar  topography  values  contained  in  the  higher- level  data  product  “CLEM1-L- 

LIDAR-5-TOPO-Vl.O”  and  to  the  Clementine  Topography  Map  of  the  Moon  (Figure  6) 
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to  determine  erroneous  laser  measurements.  The  resulting  lunar  point  cloud  of  valid 
topography  values  contains  approximately  73,000  points. 

Figure  107  superimposes  the  Clementine-derived  2D  point  cloud,  exported  from 
Cloud  Compare,  onto  a  cylindrical  lunar  map  with  a  resolution  of  64  pixels/degree.  The 
underlying  lunar  mosaic  is  one  of  the  wide-angle-camera  global  morphological  maps 
included  in  the  “Lunar  Reconnaissance  Orbiter  Camera  Reduced  Data  Record,”  courtesy 
of  NASA/GSFC/ Arizona  State  University.  The  image  is  centered  at  zero  degrees  north 
latitude  and  180  degrees  east  longitude,  corresponding  to  the  farside  of  the  Moon. 
Topography  values  are  in  meters  above  a  lunar  spheroid  with  an  equatorial  radius  of 
1,738  km. 


Figure  107.  Clementine  LiDAR  Topographic  Point  Cloud.  Adapted  from 
Berwick  et  al.  (2007)  and  Robinson  (2010). 


Next,  we  convert  the  Clementine  point  cloud  from  IAU  2000  longitude/latitude 
coordinates  to  Cartesian  coordinates,  in  accordance  with  Table  14,  to  produce  a  3D  lunar 
point  cloud.  Unfortunately,  the  point  total  of  the  Clementine  data  (-73,000  points)  is  far 
below  that  of  the  MOLA  data  (-590,000,000  points).  Thus,  the  resulting  Clementine  3D 
point  cloud  has  insufficient  point  density  to  convey  a  complete  spherical  representation  of 
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the  lunar  surface,  unless  the  view  screen  is  zoomed  out  substantially.  Despite  the  poor 
resolution,  the  lower-quality  3D  point  clouds  are  shown  in  Figure  108. 


Figure  108.  Clementine  LiDAR  Point  Cloud  in  3D — Various  Perspectives. 

Topography  Farside  (Close-up  view) 


Nearside  (Distant  view)  Farside  (Distant  view) 


From  Figure  107  and  the  “farside  (close-up  view)”  in  Figure  108,  we  see  that 
Clementine’s  approximately  350  mapping  orbits  did  not  provide  complete  coverage  of 
the  Moon’s  surface.  Future  work  could  apply  this  process  to  Lunar  Reconnaissance 
Orbiter  (LRO)  data  collected  by  its  LiDAR,  the  Lunar  Orbiter  Laser  Altimeter  (LOLA). 
LOLA  2D  and  3D  point  clouds  can  then  be  compared  to  those  from  Clementine.  LOLA 
datasets  are  also  available  from  the  NASA  Planetary  Data  System. 
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